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Abstract 



Important applications in science and engineering, such as modeling traffic flow, seismic waves, electromagnetics, and 
the simulation of mechanical stresses in materials, require the high-fidelity numerical solution of hyperbolic partial 
differential equations (PDEs) in space and time variables. Many interesting physical problems involve nonlinear and 
anisotropic behavior, and the PDEs modeling them exhibit discontinuities in their solutions. Spacetime discontinuous 
Galerkin (SDG) finite element methods are used to solve such PDEs arising from wave propagation phenomena. 

To support an accurate and efficient solution procedure using SDG methods and to exploit the flexibility of these 
methods, we give a meshing algorithm to construct an unstructured simplicial spacetime mesh over an arbitrary sim- 
plicial space domain. Our algorithm is the first spacetime meshing algorithm suitable for efficient solution of nonlinear 
phenomena in anisotropic media using novel discontinuous Galerkin finite element methods for implicit solutions di- 
rectly in spacetime. Given a triangulated (i-dimensional Euclidean space domain M (a simplicial complex) and initial 
conditions of the underlying hyperbolic spacetime PDE, we construct an unstructured simplicial mesh of the (d+l)- 
dimensional spacetime domain M x [0,oo). Our algorithm uses a near-optimal number of spacetime elements, each 
with bounded temporal aspect ratio for any finite prefix M x [0, T] of spacetime. Unlike Delaunay meshes, the facets 
of our mesh satisfy gradient constraints that allow interleaving the construction of the mesh by adding new space- 
time elements, and computing the solution within the new elements immediately and in parallel. Our algorithm is an 
advancing front procedure that constructs the spacetime mesh incrementally, in a fashion similar to the Tent Pitcher 
algorithm of Ungor and Sheffer dUSOOl l. later extended by Erickson et al. (IEGSU02I) . We extend previous work which 
was limited by conservative assumptions; in doing so, we solve the algorithmic problem of predicting discontinuous 
changes in the geometric constraints that must be satisfied by the mesh. 

In 2D X Time, our algorithm simultaneously adapts the size and shape of spacetime tetrahedra to a spacetime error 
indicator. We build on and extend the mesh adaptivity technique (lACE+04b of Abedi and others including this author. 
We are able to incorporate more general front modification operations, such as edge flips and limited mesh smoothing. 
Our new algorithm is the first to simultaneously adapt to nonlinear and anisotropic response of the underlying medium, 
thus unifying mesh adaptivity with earlier work by the author. The key insight of our unified algorithm is to exploit 
the geometric structure of the so-called cone constraints that impose gradient constraints on certain facets of the mesh. 
We present our algorithm as a two-parameter optimization problem where the objective is to maximize the height 
of spacetime tetrahedra. The parameters that guide the algorithm can be chosen adaptively by the algorithm at each 
step. Our unified algorithm can choose to refine the mesh to increase the temporal aspect ratio of spacetime elements, 
independently of the numerical error indicator Our new adaptive algorithm can better optimize the shape of spacetime 
elements and is therefore promising to produce better quality meshes. 

The front is coarsened by merging adjacent facets, whether or not they were created by a previous refinement. We 
show how to make tent pitching less greedy about maximizing the progress at each step so that a subset of the front 



can be made coplanar for coarsening in the next step. This algorithm can be used in arbitrary dimensions to make the 
entire front conform to a target time T > 0. The main contribution of this extension is to retain the progress guarantee 
of earlier algorithms up to a constant factor. 

In higher spatial dimensions t/ > 3, we are able to guarantee limited mesh smoothing to improve in practice the 
quality of the mesh. Though we are unable to extend our adaptive algorithm to spatial dimensions d > 3, we state the 
major open problem that remains. 

Empirical evidence confirms that improving the quality of the spatial projection of each front leads to better quality 
spacetime elements. We give an algorithm in arbitrary dimensions to smooth the front as by pitching inclined tentpoles, 
and by performing edge flips and edge contractions in 2D x Time, so that the quality of future spacetime elements can 
be improved in practice. We strengthen the progress constraints to anticipate future changes in the front geometry. This 
algorithm represents recent progress towards a meshing algorithm in 2D x Time to track moving domain boundaries 
and other singular surfaces such as shock fronts. 

Additionally, we are able to prove theoretical bounds on the worst-case temporal aspect ratio of spacetime ele- 
ments, while at the same time providing parameters to the algorithm that improve either the number or the aspect ratio 
of elements in practice. The number of elements in our mesh is provably close to that in a size optimal mesh that also 
satisfies certain constraints required of any mesh suitable for solution by SDG methods. 
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Chapter 1 

Introduction 



In this thesis, we study spacetime meshing to support fast and accurate simulation of certain physical phenomena by 
Spacetime Discontinuous Galerkin (SDG) finite element methods. We give algorithms to construct spacetime meshes. 
We prove guarantees on the performance of the algorithm. We prove worst-case bounds on the size and quality of the 
individual elements in our meshes. We want our meshes to be efficient and lead to accurate solutions in practice. Our 
meshes support an efficient solution strategy by SDG methods because they satisfy certain geometric constraints. Our 
meshing algorithms exploit the flexibility of novel SDG methods. 

Mesh generation is an important problem in computational geometry. Delaunay meshes jdBSvKOOOI l are widely 
used in engineering applications. Solving engineering problems using a mesh imposes geometric requirements on 
the shape and size of mesh elements (IShe02l l. The efficiency of the solution technique depends on the number and 
distribution of elements in the mesh. Delaunay meshes are known to nicely satisfy many such requirements in practice 
and have very good theoretical properties such as size optimality. Our spacetime meshes are not Delaunay meshes but 
our spacetime meshing problem is motivated by similar geometric requirements and by a similar desire for efficiency. 
We confront and solve many challenges unique to meshing a non-Euclidean spacetime domain. We do so by building 
on previous research JUSOOi |U02[ IEGSU021 lEGSU05l l and insight due to several researchers. The research in this 
thesis is part of a larger project at the Center for Process Simulation and Design (CPSD) at the University of Illinois 
at Urbana-Champaign. 

The algorithms in this thesis are being implemented by members of CPSD, including the author, and tested with 
actual DG solvers on practical simulation problems to illustrate the usefulness and applicability of the research. The 
reader is referred to other publications ( iACE+04[ lACF+041 IAHTE051 ) for empirical results. 

1.1 Background 

In this section, we will review basic concepts about hyperbolic partial differential equations (PDEs) and physical 
phenomena modeled by hyperbolic PDEs, to motivate the application of this research to solving simulation problems 
in science and engineering. The reader is referred to classical textbooks (IStr86l IHea02t |Fol95t IWhi74l) for a review 
of the theory of PDEs, hyperbolic problems, and the finite element method. There are also several online resources. 
Jim Herod of Georgia Tech has notes and Maple worksheets online (IHer04l) . See also Math World (IWei99bl) and links 
therein. Introductory texts on PDEs include those by Tveito and Winther (ITW98I I. Folland (IFol95l) . Cooper (ICoo98l l. 
and Thomas ( |Tho99l l. 

Simulation problems in mechanics consider the behavior of an object or region of space over time. Scientists 
and engineers use conservation laws and hyperbolic partial differential equations (PDEs) to model transient, wave- 
like phenomena propagating over time through the domain of interest. PDEs are mathematical models of physical 



phenomena, often derived by applying a fundamental principle such as conservation of mass, momentum, or energy. 
Example applications are numerous, including, for instance, the equations of elastodynamics in seismic analysis and 
the Euler equations for compressible gas dynamics. Computational mechanics is concerned with the accurate computer 
simulation of physical phenomena. 

An example of a hyperbolic PDE is the one-dimensional wave equation 

2, 



d u r,d u 



= o?^-^ forO<x<Landf >0 



(9f2 dx^ 



with boundary conditions 



and initial conditions 
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du 

Tt 



r=0 



This PDE describes, for instance, the vibration of a taut string when it is plucked. Here, u{x,t) is the displacement of 
point X at time t and co, a constant, is the wavespeed, i.e., the speed at which the wave travels along the string. 

A PDE is homogeneous if the trivial function m = is a solution. 

A PDE is linear if the PDE and the initial or boundary conditions do not include any product of the dependent 
variable u or its partial derivatives. A PDE that is not linear is nonlinear. For example, 

• First-order linear PDE: m, + CU;^ = 

• Second-order linear PDE: Uxx + Uyy ~ (j) {x,y) for an arbitrary function (p of the independent variables only. 

A linear PDE has the property that if mi and U2 both satisfy the PDE, then so does every linear combination aui+j5u2 
This property is called the principle of superposition. 

A nonlinear equation is semi-linear if the coefficients of the highest-order derivatives are functions of the indepen- 
dent variables xi, X2, . . ., x„ only. Examples of semi-linear PDEs are the following: 



{x + 3)ux+xy Uy = u 



3 



• 



xuxx + {xy + y ) Uyy + uux + u Uy ^ u 



A nonlinear PDE of order m that is not semi-linear is quasi-linear if the coefficients of the derivatives of order m 
depend only on the independent variables and on derivatives of order less than m, e.g., 

(1 + Uy)Uxx - lUxUyUxy + (1 + u\)Uyy = 

A nonlinear PDE of order m is fully nonlinear if the highest-order derivatives of u appear nonlinearly in the 
equation. 

The following examples illustrate the distinction between the different types of PDEs: 

• Mr + Mjc = is homogeneous linear 

• Mr + xux =: is homogeneous linear 

• Uxx -\-Uyy = Q is homogeneous linear 



• Mf+jc^M;^ = is homogeneous linear 

• Uf + Uxxx + MM.T = is homogeneous semi-linear 

• Ut + Ux + u^ ~Q is homogeneous semi -linear 

• M„ + Uyy —x^+y^ is inhomogeneous linear 

• M, + uux = is not semi-linear, but it is quasi-linear 

• M^ + My = 1 is inhomogeneous fully nonlinear 

Many interesting physical problems are modeled by second-order PDEs. The general form of a linear second-order 
PDEis 

oUxx + buxt + cuii + dUx + eut + fu + g(x,t) =0 

This PDE is classified as (i) parabolic if h^ — Aac — 0, (ii) hyperbolic if b^ — Aac > 0, and (iii) elliptic if b^ — Aac < 0. 

1.1.1 Numerical methods 

It is very difficult and often impossible to obtain analytical solutions of PDEs arising in practical problems. Therefore, 
approximate solutions are obtained by numerical techniques. The choice of numerical method depends on the physical 
behavior of the system, i.e., on the type of PDE that models the system. 

Several numerical techniques can be used for numerical solutions of partial differential equations. They are clas- 
sified on the basis of the discretization method by which the continuum mathematical model is approximated. These 
include the finite difference method, finite volume method (FVM), meshfree methods, and finite element method 
(FEM). 

The finite element (FEM) method is a vastly popular numerical technique for solving PDEs. The FEM method 



requires a mesh of the domain. There is a vast body of literature on mesh generation. In Section [TTL2] we will briefly 
mention the results most relevant to our problem to place our meshing algorithms in context. 

Discrepancies can arise between the actual physics and the approximate numerical solution of the PDE modeling 
the physics. The sources of such discrepancies include errors in modeling physical phenomena as well as errors due 
to computation such as discretization error and numerical roundoff error 

Some errors are caused because of the discretization of the domain. Errors are also caused by truncation and 
rounding due to finite precision of the computations. The order of accuracy of a numerical method is the exponent 
of the largest order term in the Taylor series expansion of the difference between the analytical solution and the 
approximate solution. Numerical instability occurs when the error term grows so fast that it dominates the actual 
solution. Both accuracy and stability affect the computational time required to converge (if at all) to the approximate 
solution. Therefore, it is important to minimize error and maximize stability. Since there is always a tradeoff between 
accuracy and efficiency of any numerical method, the most desirable methods are those for which the running time 
increases slowly as the acceptable error decreases. 

Galerkin method 

The Galerkin method ( IWei99al l provides a unifying basis from which several numerical methods, like finite element 
and finite difference methods, can be derived. The Galerkin approximation to the analytical solution of a PDE is an 



approximate numerical solution that is the unique best approximation in a certain sense to the exact solution. We refer 
the reader to other resources, for instance the online lecture notes by the Chalmers Finite Element Center (ICen03b . for 
a detailed discussion of the Galerkin method. 

Discontinuous Galerkin method 

The discontinuous Galerkin (DG) method, originally developed by Reed and Hill in 1973 for neutron transport prob- 
lems and first analyzed by Le Saint and Raviart in 1975, is used to solve ordinary differential equations and hyperbolic, 
parabolic, and elliptic partial differential equations. 

The DG method is a technique that uses discontinuous basis functions to formulate a Galerkin approximation. 
Given a mesh of the analysis domain, the DG method approximates the solution within each element by a function 
from a low-dimensional vector space of functions, e.g., as a linear combination of basis functions like polynomials. 
For a pair of adjacent mesh elements, the approximate solution computed in the interior of the elements does not have 
to agree on their common boundary. 

The DG method has many desirable properties that have made it popular For example (i) it can sharply capture 
solution discontinuities relative to a computational mesh; (ii) it simplifies adaptation since inter-element continuity is 
neither required for mesh refinement and coarsening, nor for p-adaptivity; (iii) it conserves the appropriate physical 
quantities (e.g., mass, momentum, and energy) on an element-by-element basis; (iv) it can handle problems in complex 
geometries to high order; (v) regardless of order, it has a simple communication pattern to elements sharing a common 
face that simplifies parallel computation. If the solution field has a discontinuity in the form of a shock, SDG methods 
can exactly capture this discontinuity if the mesh facets are perfectly aligned with the shock surface. On the other hand, 
with a discontinuous basis, the DG method produces more unknowns for a given order of accuracy than traditional 
finite element or finite volume methods, which may lead to some inefficiency. The DG method is harder when applied 
to unstructured meshes; in particular, it is harder to formulate limiting strategies to reduce spurious oscillations when 
high-order methods are used. 

For a thorough discussion of DG methods and a review of the state of the art, we refer the reader to the book by 
Cockburn, Kamiadakis, and Shu JCKSOOb . 

Spacetime discontinuous Galerkin methods 

The first discontinuous Galerkin formulation for spacetime problems, introduced by Reed and Le Saint et aJ., was 
discontinuous in time only. The spacetime domain was divided into slabs by constant-time planes. The solution 
was assumed to be continuous within a slab and jump conditions at the inter-slab boundaries were used to enforce 
the appropriate level of continuity between adjacent slabs. The spacetime basis functions of this so-called time- 
discontinuous Galerkin method were continuous in spacetime except across the constant-time planes that separated 
two consecutive slabs. The spacetime mesh was, therefore, required to conform to the inter-slab boundaries. 

Richter (IRic94l l introduced the first spacetime-discontinuous Galerkin finite element method, which permitted ba- 
sis functions to be discontinuous in spacetime across all element boundaries, for the wave equation. Yin etal. ( |YAS+99l 
lYAS^OOl IYin02b proposed the first spacetime-discontinuous Galerkin method for linear elastodynamics. Abedi et 
al. ( IACF+04[ IAHP051 lAHTEOSt give an improved formulation. Since the SDG method allows every pair of adjacent 
elements to have basis functions that are discontinuous across their common boundary, the SDG method can solve 
fully unstructured and even nonconforming meshes. The flexibility of SDG methods, in terms of the types of mesh 
that are solvable, is exploited by the meshing algorithms in this thesis. 



1.1.2 Mesh generation 

Given a domain of interest, the mesh generation problem is to decompose the domain into simple cells, called elements, 
each of constant complexity. Such a decomposition is called a mesh. Elements consist of vertices, edges, ridges, and 
facets of appropriate dimensions. Linear, bilinear, or trilinear elements are often used; convex cells (e.g., simplices) 
are sometimes preferred. A triangulation is a decomposition into simplices ( |Zie95l l. In general, a mesh may contain 
hexahedra, prisms, pyramids, and other simple shapes, possibly in the same (hybrid) mesh. 

If the domain is bounded, the mesh must conform to the boundary, i.e., the domain boundary must be contained in 
the union of a subset of the mesh faces. A curved domain boundary cannot be captured exactly by a mesh consisting 
of linear elements; in this case, the mesh must conform to the boundary approximately in a piecewise linear fashion. 

An important question is whether a domain can be meshed at all and, if so, with how many elements. It is well- 
known that every simple polygon in the plane can be triangulated (IGJPT78 i in fact in linear time ( ICha9 1 1 lAGROOI l . In 
dimensions three and higher, additional vertices, called Steiner points, may be required; in general, it is NP-complete 
to determine if one or more Steiner points are required and to compute the minimum number of Steiner points ( |RS92| l. 

The accuracy and stability of numerical methods are affected by the geometric shape of mesh elements. Therefore, 
it is important to generate a mesh with only good quality elements. The notion of element quality depends on the 
numerical method (IShe021 l. 

When the domain is Euclidean, several popular metrics are used to define the quality of the mesh elements and of 
the mesh as a whole. Popular objectives are to optimize at least one of the following quality measures for every element 
in the mesh: minimize ratio of circumradius to length of shortest edge, maximize smallest angle, minimize largest 
angle, etc. Different metrics appear to work better than others depending on the application. See Shewchuk (IShe02l) 
for a comprehensive study of quality metrics. In 2D, all these metrics are equivalent to each other up to constant factors; 
this is not true in higher dimensions. In 2D, Delaunay triangulations JdBSvKOOOI l have guaranteed quality JCDE+OOl 
ICD03]) . 

For anisotropic problems, quality metrics should also be anisotropic. Typically, the local anisotropic nature of the 
problem is used to define a local quality measure at each point in the domain. The problem now is to generate a mesh 
such that for every point in the domain, the quality of the element that contains the point is at least as good as the 
desired quality metric at that point. Labelle and Shewchuk ( ILS03I ) propose a definition of an anisotropic Delaunay 
triangulation and an algorithm to construct such a triangulation of guaranteed quality. 

Bern, Eppstein, Gilbert ( IBEG94 ) use quadtrees to construct triangulations with extra Steiner vertices added to the 
original input. They present the first algorithm to triangulate a planar point set and polygons, with all angles bounded 
away from zero, using a number of triangles within a constant of optimal. In higher dimensions, they are able to 
triangulate point sets such that the resulting triangulation has no small solid angles and the number of its simplices 
is within a constant factor of the optimal number The result of Bern et al. was generalized to higher dimensions by 
Mitchell and Vavasis dMVOOI l. 

A (i-dimensional structured mesh, also known as a (structured) grid, is a cubical complex (lZie95l l with the property 
that every A:-dimensional cell belongs to exactly 2'^^^ li-dimensional cells; only cells on the boundary of the complex 
are excepted from this rule. Thus, all interior nodes are incident on the same number of elements. For example, a 
structured quadrilateral grid in 2D has each interior node incident on four quads; a structured hexahedral grid in 3D 
has each interior node incident on eight hexahedra. Because the combinatorial structure of such a grid is so regular, 
simple data structures suffice to store and manipulate the grid. Unstructured meshes relax the uniformity requirement 
on the node degree. 



Mesh smoothing and mesh refinement are commonly used to improve the quality of a given mesh, for instance, a 
Delaunay triangulation of a given domain. 

Mesh smoothing, also called r-refinement, adjusts the locations of the vertices of the mesh while keeping the com- 
binatorics of the mesh unchanged. Smoothing serves to improve the quality of the elements at least locally (jABEg^l 
but global convergence to some optimum quality is not guaranteed. Mesh smoothing can be appHed to both structured 
and unstructured meshes. 

Mesh refinement, also called /i-refinement, is the operation of introducing additional Steiner points to alter the mesh 
and improve its quality. Delaunay refinement (IChe89t Rup95 i, is used when the objective it to obtain a Delaunay 



triangulation JdBSvKOOO ) of the domain where each element has some guaranteed quality. Delaunay refinement 
proceeds by destroying bad quality elements by inserting their circumcenters and updating the Delaunay triangulation. 
The refinement process eventually terminates but in 3D does not eliminate so-called slivers. Weighted Delaunay 
refinement ( ICD03b attempts to remove slivers in the next step. If the mesh is required to conform to a piecewise linear 
boundary and if the boundary contains an acute dihedral angle, it is not known whether weighted Delaunay refinement 
can produce a good-quality triangulation without slivers. 

Delaunay refinement is an incremental algorithm to improve the quality of a simplicial mesh, for instance by 
inserting circumcenters of bad-quality simplices, and updating the mesh to maintain the Delaunay property. See the 
seminal papers by Chew ( |Che89l l and by Ruppert ( Rup95 1 on the quality and optimality of the triangulations resulting 



from Delaunay refinement by circumcenter insertion. For Delaunay refinement using alternatives to circumcenter 
insertion, see Rivara's longest edge bisection ( IRiv97l l. sink insertion by Edelsbrunner and Guoy (lEGOll l. and Ungor's 



Off-Center insertion algorithm (Ung04 1. 



Another popular technique to mesh a domain is to begin with a surface mesh of the boundary of the domain and 
extend it incrementally into the interior of the domain to get a volume mesh. Such an approach is called an advancing 
front method. Advancing front techniques are popular due to their simplicity, but suffer from a lack of provable 
guarantees; for instance, an advancing front algorithm may create degenerate or inverted elements ( ISev97l ), it may 
create non-triangulable volumes and fail ( |fPCS99l l, or it may require resorting to ad-hoc techniques and heuristics 
where different parts of the front collide. 

Meshes of complicated domains frequently consist of millions or billions of elements. Representing and storing 
large meshes in memory is complicated because memory in real-world computers is organized in a hierarchy ranging 
from fast online storage to slow offline devices, and the amount of fastest memory is limited. Exploiting the spatial 
and temporal locaUty exhibited by most algorithm to compactly store a large mesh in a way that supports efficient 
access is an active and exciting area of research (ICRMS03IIYLPM05t . Minimizing memory usage is of concern to us 
too. Our algorithms minimize memory requirements by storing in main memory only a small subset of the mesh, i.e., 
the elements adjacent to the advancing front, instead of the entire spacetime mesh. 

Multiprocessor computers and computer networks are available to solve large-scale problems. Parallel algorithms 
for mesh generation ( |STU04| ) and parallel numerical methods take advantage of such computer systems by sophisti- 
cated techniques. Our meshing algorithms are highly parallelizable and the discontinuous nature of the SDG method 
makes it possible to compute the solution in unrelated subdomains simultaneously in parallel. 

Automatic generation of good quality finite element meshes — with at most minimal user interaction — ^remains 
a topic of active research interest ( |Owe99l ). In the spacetime domain, we have only begun addressing the many 
challenges with many exciting avenues to explore in the future. 



1.1.3 Causality 

The phenomenon that characterizes hyperbolic problems is causality, the fairly intuitive notion of cause and effect, 
stimulus and response, commonly observed in daily life. When a pebble is dropped into a pond, the disturbance 
radiates outwards in expanding circles from the point of origin. The result is waves traveling on the surface of the 
water, traveling with a bounded wavespeed. An example in one-dimensional space is the vibration of a taut string 
fixed at both ends when it is plucked. Let L be the length of the string. The displacement u{x,t) of the point of 
the string at coordinate x about its mean position at time t is described by a second-order hyperbolic PDE, the wave 
equation m„ — CO^u^x ~ with the initial conditions m(x,0) — for every x £ [0,L] and u{0,t) — u{L,t) = for every 
f > 0. The wavespeed CO is the rate at which the initial stimulus travels to other points of the string. The propagation of 
influence with a bounded wavespeed within a spatial domain, a subset of E'', can be visualized as a relation between 
points in the spacetime domain one dimension higher, i.e., E'' x M. 

Influence and dependence 

Points in spacetime are partially ordered by causality — a point P influences another point Q, written as P -< Q if and 
only if changing the physical parameters at P could possibly change the result at Q. A point Q depends on P, written 
as Q y P, if and only if P influences Q, i.e., P ^ Q-i^ Qy P- The domain of influence of P is the set of points Q such 
that P < Q. Symmetrically, the domain of dependence of a point P is the set of points R such that P y R. We say that 
two distinct points P and Q are independent if and only if neither P ^ Q nor Q <P. 

A characteristic of a PDE through a given point x is the iso-surface along which u satisfies an ordinary differential 
equation (ODE). Thus, characteristics are curves or surfaces or hypersurfaces, of dimension n — 1 or less, in the space 
spanned by the n independent variables. Each characteristic corresponds to a characteristic equation which is an 
ordinary differential equation (ODE) or a system of ODEs. The slope of a characteristic at P is the slope of the tangent 
to the characteristic at P. The method of characteristics for solving a given PDE changes coordinates from {x^t) to a 
new coordinate system (x(5),s) in which the PDE becomes an ordinary differential equation (ODE) in the independent 
variable s. The variable s is a parameter that varies along the characteristic hypersurface. 

For a hyperbolic spacetime PDE, a characteristic hypersurface through a point P represents the flow of information 
through the space domain with time. For hyperbolic PDEs such as the wave equation m„ — aP'Uxx = the solution at 
p — {xp,tp) influences the solution at every future point in the wedge bounded by the characteristic lines through p 
with slope ±0). The information at p propagates into the future with a uniform bounded speed CO in every direction. 

A curve F = {x,y{x)) is characteristic for the second-order PDE 

auja + buxy + cuyy + dux + euy + fu + g{x,y) — 

if 

dy _b± Vb^ - 4ac 

dx 2a 

along F. The three different classes of second-order PDEs have three different types of characteristics. Hyperbolic 
PDEs have fe^ — 4-ac > 0, so Equation ] 1.1 [ has two real solutions and there are two characteristic curves through every 
point {x,y). For example, the wave equation Uxx — Uyy = is hyperbolic and has two characteristic lines given by 
y= -iix. 

The solution to a PDE can have discontinuities. If two characteristics intersect at a point x, the two values for the 



Figure 1.1: Points A and B are independent because B lies outside the cone of influence of A and A lies outside the 
cone of dependence of B 

solution M at X obtained by back-tracing the characteristics can be different. Such discontinuities are called shocks. 
Shocks can occur even when the initial conditions are smooth. 

For spacetime hyperbolic PDEs, the wavespeed everywhere is finite, i.e., information propagates through the space 
domain with a bounded wavespeed. 

Cones of influence and dependence 

The slope g{P) at a point P is a quantity computed by the numerical solver to bound from above the speed at which 
any change in the parameters at P propagates through the spacetime domain. For instance, g{P) can be computed as 
the maximum slope of every characteristic through P. The wavespeed at P, denoted by Co{P), is the reciprocal of the 
slope, i.e., 

^ ' <y{p) 

Thus, the slope (or wavespeed) at every point in the domain defines a scalar field over the domain. 

The cone of influence of P, denoted by cone+(f ), is a circular cone with apex at P, axis in the positive time 
direction, and slope equal to o{P). The cone cone+(P) is a closed full-dimensional subset of the spacetime domain. 
Symmetrically, the cone of dependence of P, denoted by cone^(P), is a circular cone with apex at P, axis in the 
negative time direction, and slope equal to o{P). See Figure [lT] 

When the underlying medium is anisotropic, due to nonlinear response, waves propagate faster in some spatial 
directions than in others. In such cases, cones of influence need not have circular cross-sections. We will postpone a 
discussion of this anisotropic situation until later in this section. 

No-focusing assumption 

We make the following assumption about how well the cone of influence cone+ (P) at a point P approximates the 
domain of influence of P. We assume that the slope of the cone of influence at a point P in the future can be estimated 
from the slopes of the cones of influence of all points Q such that P G cone+Q. In other words, we assume that P is 
influenced only by points Q such that P e cone+2- This property is formally stated as Axiom [Ll] 




Figure 1.2: No-focusing means that the wavespeed at a point in the future can be estimated from the wavespeeds of 
points on the current front. 

Axiom 1.1 (No focusing). For every point P in the spacetime domain, the slope o{P) is bounded by the minimum and 
maximum slope of every cone of influence containing P. Thus, 

min {a{Q)} < a{P) < max {aiQ)} 

P e cone+ (Q) Pe cone+ (g) 



In Figure 1.2 point P is influenced by point Q, hence C7(P) — O'(C) is possible. Suppose P is advanced in time 



to P'. The point P' is influence by both Q and Q\ Hence, g{P) — mm{a{Q),G{Q')} is possible. 

If a point P = (p, f ) is in a cone of influence C, then every point P' = {p,t') where f' > f is also in C. See Figure [L2| 
Together with Axiom [TT[ we obtain the following monotonicity lemma. 

Lemma 1.2. For an arbitrary flxed point p G E , let P — {p,t) and P' — {p,t') be two points in spacetime with the 
same spatial projection p and time coordinates t and t' respectively. Ift' > t, then cy{P') < <7(P). 

Explanation of the no-focusing assumption: Axiom [TTT] The causal slope a{P) is a first-order approx- 
imation to the slope of the characteristic curves/surfaces through P; the slope of the cone of influence at P in every 
direction is less than or equal to the slope of all the characteristics through P in that direction. Axiom [TTT] which holds 
in the absence of focusing, allows us to conservatively estimate the causal slope at every point in spacetime where the 
solution has not been computed by our algorithm yet. Cones of influence are locally conservative approximations to 
the domains of influence. A sufficiently conservative cone can be chosen that is a valid approximation of the actual 
domain of influence for the finite step size chosen by our algorithm to advance the solution. Thus, in Figure |1.1[ 
points A and B are independent because neither influences the other. 

We assume only that the no-focusing assumption of AxiomfTTTIholds everywhere in the spacetime domain £1. Our 
algorithms guarantee correctness of the solution only in the absence of focusing. 



Asymmetric cones Due to nonlinear response and anisotropy of the underlying medium, waves can propagate 
with different speeds in different directions. Therefore, in general, the wavespeed at a point P is a function of the 
spatial direction n, and hence, cones of influence and dependence have non-circular cross-sections. Our assumption 
of no-focusing in this anisotropic context can be stated as the following axiom. 

Axiom 1.3 (Anisotropic no-focusing). For every point P in the spacetime domain, the slope at P in an arbitrary spatial 
direction n, denoted by 0„(P), is bounded by the minimum and maximum slope in the same spatial direction n of every 
cone of influence containing P. Thus, for every spatial direction n, we have 

min {a„(e)} < (7„(P) < max {(7„(e)} 

Pecune+{Q) Pecone+(Q) 



1.2 Spacetime meshing 

We saw in Section [LT] that wave propagation is modeled by hyperbolic partial differential equations (PDEs) in both 
space and time variables, for instance, the wave equation m„ — (O^u^x = in ID space x time. The wavespeed (0, 
the speed at which changes in physical parameters at a point (x,f) propagate to other points in the domain, may 
be a function of x and t as well as of u and its derivatives. The spacetime discontinuous Galerkin (SDG) method 
approximates the solution within each spacetime element of a mesh of the spacetime domain as a linear combination 
of simple basis functions such as polynomials. The SDG method allows basis functions to be discontinuous across 
element boundaries — adjacent elements need not agree on the solution along their common boundary. 

The spacetime DG method motivates our meshing problem. In this chapter, we will give an overview of the goal 
of this thesis — to give provably correct algorithms for generating efficient spacetime meshes suitable for DG solvers. 
We will describe the challenges that are encountered. We will enumerate the assumptions under which our algorithms 
and the meshes they produce satisfy the stated goal. We conclude this chapter with a summary of our results and of 
previous research on spacetime meshing. 

This section defines the basic notions and vocabulary used throughout the rest of the thesis. We will remind the 
reader of our assumptions and refer back to this section whenever necessary. 

Patches 

The notions of influence and dependence extend naturally to arbitrary subsets of spacetime. For arbitrary elements A 
and B, not necessarily distinct, of a spacetime mesh, we say that A influences B, written as A ^ B if and only if some 
point a E A influences some point b E B. If A ^ B, then A must be solved no later than B, otherwise the solution in 
B is not valid. Elements A and B are coupled if and only if both A ^ B and B ^A. A pair of coupled elements must 
be solved together, i.e., as a simultaneous system of equations. If neither A ^ B nor B ^ A, then we say A and B are 
independent. 

A patch n in a spacetime mesh i2 is a set of elements that must be solved simultaneously. More formally, a patch 11 
is a set of one or more elements with the following properties: 

1 . for every A G 11 there exists B e 11 such that A and B are coupled (dependency condition); 

2. for every A G 11, if there exists B G ii such that A and B are coupled then B G 11 (maximality condition). 

10 



The size of a patch is the number of elements in it. 

We have seen that causaHty defines a partial order over points in spacetime. We say that a subset Q! of spacetime is 
causal if and only if Of is an anti-chain in the partial order, i.e., if and only if no two points of Df influence each other 
The boundary of a patch n is necessarily causal because otherwise elements outside the patch would be coupled with 
elements in the patch, violating the maximality of H. Patches are therefore partially ordered by causal dependence. 

A spacetime mesh Q. can be solved patch-by-patch in an order that respects the partial order of patches. The 
total computation time of such a patch-wise solution strategy is the sum over every patch n in the mesh Q. of the 
time to solve Yl. We say that a patch-wise solution strategy is efficient if the total computation time is proportional 
to the number of elements in the mesh. Our meshing algorithms support an efficient solution strategy because both 
(i) the maximum size of a patch and (ii) the number of patches in a given spacetime volume are bounded. In a mesh 
constructed by other algorithms, either the maximum size of a patch or the number of patches or both are unbounded; 
furthermore, the solver must explicitly compute the partial order of patches or solve all patches together as a large 
coupled system. Therefore, a generic purpose meshing algorithm is inefficient in general; specialized algorithms are 
necessary for efficiently meshing in spacetime. 

1.2.1 Terminology and notation 

Before we proceed to the technical details of our spacetime meshing problem and a summary of the results in this 
thesis, we define some basic concepts and introduce the reader to the notation used in the rest of this thesis. 

First, we review some basic concepts. See the references (IZie95b for further details. 

A simplex of dimension k, also called a fc-simplex, is the convex hull of A:+ 1 affinely independent points. A 
simplex is a closed set of points. We identify a simplex with its vertices. Let popip2- ■ -Pk be a A:-simplex. Every 
subset of {/7o,pi,P2i- --jPyt} defines a /flce of the simplex. Simplices of dimension 0, 1,2, and 3 are called vertices, 
segments, triangles, and tetrahedra respectively. For a A:-simplex S, let affS' denote the affine hull of 5, i.e., the set of 
all linear combinations of vertices of 5; the affine hull aSS is a (^— \)-flat JdBSvKOOOl l. For example, the affine hull 
aifpq of two points p and q is the line pq. 

A simplicial complex is a collection ^^ of simplices with the following properties: (i) if a simplex S £'^ then 
every face of S is in '^, and (ii) two simplices in "^ either do not intersect, or their intersection is a simplex of smaller 
dimension which is their common face of maximal dimension. The empty simplex, whose dimension is — 1, is a face of 
every simplex. The dimension of the simplicial complex "^ is the highest dimension of any simplex in the collection '^. 
For a (i-dimensional simplicial complex, the {d~ 1) -dimensional faces are called i\\t facets. A simplicial complex is 
a special type of cell complex ( IZie95l) . 

Two faces of a complex are incident if one is included in the other Given a simplicial complex ^, the star of a 
face F, denoted by star(F), is the sub-complex consisting of aU simplices incident on F, and all their faces. The star 
of F is a closed set. The link of a face F, denoted by link(F), is the sub-complex consisting of all faces G of simplices 
in star(F) such that GHF = 0. 

A particular type of simplicial complex is a triangulation. A d-dimensional triangulation is a ^/-dimensional 
simplicial complex such that all maximal faces are ^/-simplices. Note that a triangulation need not form a manifold. 
In a J-triangulation, every {d— l)-face F belongs to one or more li-simplices; in the former case, F is a boundary 
facet, and in the latter case, F is an interior facet. All faces of a boundary facet are boundary faces; all other faces 
are interior faces. For example, if v is a vertex of a 2-dimensional triangulation, then star(v) consists of the vertex v, 
all edges incident on v, and all triangles incident on v together with their edges and vertices; link(v) is a set of edges 
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which form a closed loop if star(v') is homeomorphic to a disk. 

Next, we define the terms and notation used to describe our spacetime meshing problem. 

A (i-dimensional space mesh is a ^-dimensional triangulation. We visualize the {d+ 1) -dimensional spacetime 
domain with the time axis drawn vertically so that time increases upwards. We use uppercase letters like P, Q, R to 
denote points in spacetime and corresponding lowercase letters like p, q, r to denote their spatial projections. 

h front is a maximal subset of spacetime that forms an anti -chain in the partial order of dependence; i.e., a front T 
is a maximal set of points such that no two points of T influence each other An example of a front is the set of all 
points with a given time coordinate T, whose graph is the horizontal constant-time plane t = T. 

The front X defines a piecewise hnear function T : E'^' ^ ]R. Equivalently, the front T is a (i-dimensional piecewise 
linear terrain JdBSvKOOOI l. a subset of E"' x M. We will not distinguish between these two equivalent descriptions of 
the front — one as a function and the other as a set of points. Each point P on the front T can be written as P = (p, T{p)) 
where p is the spatial projection of P. For a real T and front T, we use T > T to mean that every point P = {p, T{p)) 
of T satisfies 'i{p) > T, i.e., the entire front T has achieved or exceeded time T. Arbitrary fronts T and t' satisfy t' > T 
if and only if T'{p) > T{p) for every point p such that P = {p, T{p)) e T and P' = {p, i'{p)) E t' ■ 

We will frequently assume that the front T consists of a single ^-simplex. In this case, we will not distinguish 
between the simplex and the corresponding time function T : E*^ ^ M whose graph is the fe-flat in spacetime spanned 
by the simplex. In such a scenario, since T is a linear function, its gradient, denoted by Vt, is the same everywhere — 
V T is a vector in M*^ in the direction of steepest ascent and its L2-norm is denoted by || V t|| . 

The front at the beginning of each iteration of our algorithm represents the frontier of both the incremental mesh 
construction as well as the solution. Thus, if P= {p,T{p)) is apoint of the front T then p belongs to the space domain 
at time T{p), i.e., p EM^-fpy, the solution at f has been computed; and the set of points f+ = {p,t^) such that f+ > T{p) 
remain to be meshed and the solution there is unknown. 

A front is causal if and only if for every facet F of the front all characteristics cross F from one side to the other 
For an element with a facet on the front T, this facet is an outflow facet of the element because influence travels 
from the interior of the element to the exterior across this facet. Equivalently, a front T is causal if and only if every 
point P = {p, t{p)) depends only on points below the front, i.e., if P y Q then Q can be written as Q= {q,tg) where 
tci < l{q). Symmetrically, every point of the front T influences only points above the front. By our earlier definition, 
every front is automatically causal. 

For a simplex (of any dimension) S of the mesh M, let t|^ denote the time function T restricted to S and extended 
to the affine hull of S; in other words, TJ^ is a linear function that coincides with T for every point of S. 

Let T, : M ^ M denote the front after the /th step of the algorithm; To is the initial front. At every step T, -^ T,+i, 
our algorithm will exphcitly triangulate the new front T,+i. Thus, every front constructed by our algorithm will be a 
triangulated terrain. For every /, the front T, is a terrain whose facets are li-simplices. 

A local minimum of the front T is a vertex p such that T{p) < T{q) for every vertex q that is a neighbor of p. 
Every front has a local minimum because it has at least one global minimum which is also a local minimum. When 
the current front T is clear from the context, for every point p e M we use P to denote the corresponding point on the 
front, i.e.,P= {p,T{p)). 

For an arbitrary simplex PqP\P2. . .Pd of the front T, we say that a vertex, say Pq, is a lowest vertex of the simplex if 
t(po) < niino<,<f/ l{pi)- Note that a simplex may have one or more lowest vertices. We say that a facet, say PiP2. ■ -Pd 
is a highest facet if the opposite vertex Pq is a lowest vertex. Note that a simplex may have more than one highest 
facets. 
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Figure 1.3: The temporal aspect ratio of a spacetime element is the ratio of its height to its duration. The height of 
both spacetime triangles in the figure is equal to the height of the tentpole they share. The spacetime triangle on the 
right has a better (larger) temporal aspect ratio than the triangle on the left. 

We say that a front t' is obtained by advancing a vertex P = {p, T{p)) of T by A/ > if T'{p) = x{p) + Af and for 
every other vertex q^ pv^s have T'{q) — T{q). For every front T, vertex p, and real A/ > 0, let t' = advance(T,/?, Af) 
denote the front obtained from T by advancing p by Af; the functions T and t' coincide everywhere outside star(/:i). 

For a simplex A, we use (t(A) to denote the minimum slope o{P) over all points P of A; we let ©(A) denote 
l/(j(A). Let (Jmin denote min/>gn{(j(P)} and ffmax denote maxpgn{ff(f )}. We assume that < (Jmin < O'max < °°- 
Thus, the slope (y{P) at an arbitrary point P is bounded: < Omm <<^{P) < O'max < °°- 

For a vertex P of the front T, define the quantity Wp in the spatial projection as 



Wn'-— min dist(D,affA). 

Aelinkip) 



For an arbitrary subset K of E'' x M, the height of K is the length of the longest time interval contained in the 
closure ofK, and the duration of /T is the length of the shortest time interval containing the interior ofK. The temporal 



aspect ratio of a spacetime element is the ratio of the height of the element to its duration. See Figure 1.3 Note that 
the temporal aspect ratio is always in the range (0, 1] with a larger value corresponding to a "better" element. The 
temporal aspect ratio of an arbitrary subset of spacetime is unchanged if all time coordinates are uniformly translated 
and scaled by a non-zero constant. 



1.2.2 Meshing objectives 

We want a mesh of the spacetime domain Q. consisting of spacetime simplices or elements such that every spacetime 
element has at least one causal outflow facet, a necessary condition for the numerical problem to be well-posed (IJJ04I I. 
The sizes of the spacetime elements are controlled by a spacetime error indicator. The spacetime elements must form 
a weak simplicial complex; this property allows a simple strategy to compute integrals over the common intersection 
of any two spacetime elements — the intersection is an entire face of at least one of the two elements and has a nat- 
ural parameterization that facilitates integration. Spacetime DG solvers are particularly suited to solve such weakly 
conforming meshes because of their discontinuous formulation. 
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Definition 1.4 (Weak simplicial complex). A weak simplicial complex is a collection 'lo of simplices with the following 
properties: 

1. if a simplex S € '^, then every face ofS is in '^; and 

2. the intersection of every two simplices o/*^ is a face of at least one of them. 

We have seen that an efficient solution strategy is to solve the mesh ii patch-by-patch in an order that respects the 
partial order of patches. Every spacetime mesh can be partitioned into patches — giving a causal partition — and can be 
solved by this strategy. The total computation cost is the sum of the cost of solving each individual patch. A generic 
mesh of spacetime cannot be efficiently solved because either the size of a patch or the number of patches or both are 
not bounded. Our objective is to minimize the total computation cost. We give a meshing algorithm that incrementally 
constructs the spacetime mesh patch-by-patch explicitly in the order of causal dependence. We prove that both the 
size of each patch and the total number of patches in the mesh of a given spacetime volume are bounded; hence, the 
total computation cost is bounded. 

We give an advancing front algorithm such that for every T G M-" there exists a finite integer k>Q such that 
the front Tj^ after the k\h iteration of the algorithm satisfies T^. > T. The input to our problem is the initial front Tq 
and the initial conditions of the PDE. We obtain the spacetime mesh by triangulating the spacetime volume, called a 
tent, between each pair of successive fronts T, and T,+i such that each spacetime element has at least one facet on the 
front T,+i, which by definition is causal. The elements in a tent are causally dependent and must be solved as a coupled 
system. We produce an efficient mesh such that the number of elements in a coupled system is small, depending on 
the maximum degree of the initial space mesh. 

Our mesh is also efficient in the sense that elements are non-degenerate in the sense that each element has temporal 
aspect ratio bounded from below. We expect that the non-degeneracy guarantee means that in practice the numerical 
solution within each element can be computed with acceptable accuracy. 

Our goal is to support an accurate solution of the underlying PDE as well as to reduce total computation time. It 
is challenging enough to devise meshing algorithms that provide a theoretical guarantee that the meshing algorithm 
would terminate after creating a finite number of non-degenerate elements. In addition, our goal is to provide an 
algorithm that is easy to implement and that performs significantly better in practice than the theoretical guarantee, for 
instance by a good choice of parameters to the algorithm and possibly using different heuristics. At the same time, we 
will never compromise on the ability to prove correctness of all algorithms described in this thesis. Some aspects of 
the meshing algorithm go beyond this basic goal and improve its performance in practice. 

Advancing front spacetime meshing 

Given a simplicial mesh of some bounded domain M C E"^, we give a meshing algorithm to incrementally construct 
a simplicial mesh of the spacetime domain using an advancing front method. The algorithm advances the front by 
moving a vertex forward in time, thus also advancing the local neighborhood star(v), and adding simplices in the 
volume between the old and the new fronts. The inflow and outflow boundaries of each patch (Figure [T3] l are causal 
by construction, i.e., each boundary facet F separates the cone of influence from the cone of dependence for every point 



on F (Figure 1.4 1. Equivalently, for every point P on F we have || VF|| < cy{P). If the outflow boundaries of a patch 
are causal, every point in the patch depends only on other points in the patch or points of inflow elements adjacent to 
the inflow boundaries of the patch. Therefore, the solution within the patch can be computed as soon as the patch is 
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Figure 1 .4: A causal triangle separates the cones of influence and dependence at every point on the face. 

created, given only the inflow data from adjacent inflow elements and initial or boundary data where appropriate. The 
elements within a patch are causally dependent on each other and must be solved as a coupled system. Patches with 
no causal relationship can be solved independently. To minimize undesirable numerical dissipation and the number of 
patches, we would like the boundary facets of each patch to be as close as possible to the causality constraint without 
violating it. 

We require that for every target time value T the algorithm will compute in a finite number of steps a mesh of the 
spacetime volume M x [0, T] and the solution everywhere in this volume. The target time T may not be known a priori 
because it depends on the evolving physics. 

Our algorithm ensures that every spacetime element in a patch has one facet on the new front, which is a causal 
outflow facet. 

Triangulating a tent Our advancing front algorithm repeatedly advances a local neighborhood A^ of the cur- 
rent front T to a corresponding neighborhood A^' of a new front t', where dN — dN' and T = t' everywhere in t\A^. 
For instance, our algorithm pitches a local minimum of the current front creating a new front which is the input to the 
next iteration of the algorithm. Pitching a vertex P of the front T means advancing P = [p, T{p)) to P' = {p, l'{p)) 
where T and t' are the old and the new fronts respectively. Pitching P means advancing the neighborhood A^ = star(P) 
to A^' = star(f') where dN = dN' = link(f ). The progress is defined to be the time difference T'{p) — T(p). 

The volume between the new front and the old front is called a tent and the edge PP' is called the tentpole. The 
tent is partitioned into spacetime elements (simplices); the set of elements in a tent forms a patch. Let F' denote an 
arbitrary facet on the new front t'. The tent is triangulated so that for each such facet F' the patch contains the element 
corresponding to the convex hull of F' and the vertex P, the bottom of the tentpole. 

For each point U in the tent, the ray PU intersects a simplex F' on the new front t'. Therefore, U is contained in 
the convex hull of F' U {P}. Hence, the set of elements in a patch partition the tent. 

By construction, each element E in the patch is the convex hull of F' U {P} for some facet F' on the new front t'. 
By construction, F' is a causal facet of E. For every point Q on F' the cone of influence of Q is entirely in the future 
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Figure 1.5: Cross-section of a patch of tetrahedra; the inflow and outflow faces are causal. Time increases up the page. 

and does not intersect the element E; therefore, F' is an outflow facet of E. Hence, every element in a patch has a 
causal outflow face. 

Also, the intersection of every two elements is either (i) empty, (ii) a common vertex, (iii) the tentpole of a patch 
containing both elements, or (iii) the entire facet of at least one of the two elements. If the two elements belong to the 
same patch with tentpole PP', then they share the vertex P; they also share either the tentpole PP' or a common facet, 
or they share a common implicit facet containing the tentpole. On the other hand, if one element is an inflow into 
the other tetrahedron then their common intersection is either the entire outflow face of the first element or the entire 
inflow face of the second element. 

When we generalize to other operations that advance the front we will ensure that these properties are always 
satisfied. For instance, an edge flip in 2D x Time creates a patch with a single tetrahedron, which has two causal inflow 
facets and two causal outflow facets. 

1.2.3 Meshing challenges 

We say that a front T is valid if there exists a positive real 5 bounded away from zero such that for every T e M-" there 
exists a sequence of fronts T, Ti , T2, . . ., T^ where T^ > T, each front in the sequence obtained from the previous front by 
advancing some vertex by 5. The minimum value of 5 for which we can guarantee that our algorithms construct only 
valid fronts is the progress guarantee of our algorithm. The progress guarantee is the same as the worst-case tentpole 
height. All our theoretical results, such as finite termination, worst-case temporal aspect ratios of spacetime elements, 
and near-size-optimality follow from the progress guarantee. What makes the definition of a valid front nontrivial is 
the requirement that all fronts be causal. The main difficulty in characterizing valid fronts arises when the wavespeed 
at a given point in the space domain increases discontinuously and unpredictably over time. In spatial dimensions 
c/ > 2, an additional complication is the changing front geometry due to adaptive refinement and coarsening, and mesh 
smoothing. 
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1.2.4 Summary of results 

The main contribution of our work is an efficient incremental spacetime meshing algorithm that adapts quickly to 
changing geometric constraints, that supports an efficient parallelizable solution strategy, and that has provable worst- 
case guarantees on the temporal aspect ratio of spacetime elements. 

Theorem 1.5. Suppose we are given a simplicial mesh M G E as well as initial and boundary conditions of the 
underlying spacetime hyperbolic PDE. Our algorithm builds a simplicial mesh Q. of the spacetime domain M x [0, °°) 
that satisfies all the following criteria: 

1. the elements of £l form a weak simplicial complex that conforms to the initial mesh M; 

2. each spacetime element has at least one causal outflow facet; 

3. for every T >0 the spacetime volume M x [0, T] is contained in the union of a finite number of simplices ofD,; 

4. the minimum temporal aspect ratio of every spacetime element is bounded from below. 

Additionally, in 2DxTit7ie, our algorithm adapts the size and duration of spacetime elements to a spacetime error 
indicator We guarantee that the diameter of each tetrahedron is no larger than that allowed by the spacetime error 
indicator Provided the error indicator does not reject any tetrahedra smaller than a bounded minimum size, our 
algorithm terminates with a finite mesh ofM x [0, T] for every target time T >0. 

Given a triangulation M of the space domain and a target time T, we say that a simpUcial spacetime mesh of 
M X [0, T] is solvable if (i) each spacetime element has both a causal inflow facet and a causal outflow facet; and 
(ii) for every point x in the spatial projection A of each spacetime element, the diameter of A does not exceed the 
diameter of the simplex of M containing x. 

Theorem 1.6. Let e be an arbitrary constant, a parameter to our algorithm, in the range < £ < j- The size of the 
mesh constructed by our algorithm is O l^j times the minimum size of any solvable mesh of the spacetime volume 
M X [0, T], where the hidden term in the 0-notation depends on the dimension d and the worst-case spatial geometry 
of each front. 



An example of a tetrahedral mesh of 2D x Time constructed by our algorithm is given in Figure [L6| Additional 
examples of meshes constructed using our algorithm appear in the references d ACE+041 lACF+041 IAHTE05 1 IThi04l l. 

1.2.5 Previous work 

This thesis is part of an ongoing long-term project at the Center for Process Simulation and Design (CPSD), a multi- 
disciplinary research group at the University of Illinois (UIUC). See |http://www.cpsd.uiuc.edu/| for more about CPSD. 
Building on ideas from earlier specialized algorithms, Ungor and Sheffer (|US02| | and Erickson et aJ. (iEGSU02l) 
developed the first algorithm to build graded spacetime meshes over arbitrary simplicially meshed spatial domains, 
called Tent Pitcher Unlike most traditional approaches, the Tent Pitcher algorithm does not impose a fixed global 
time step on the mesh, or even a local time step on small regions of the mesh. Rather, it produces a fully unstructured 
simplicial spacetime mesh, where the duration of each spacetime element depends on the local feature size and quality 
of the underlying space mesh. See Erickson et al. ( |EGSU02| | for sample meshes constructed by Tent Pitcher Figure 
from the paper by Erickson et al. (IEGSU02|) illustrates the first few tent pitching steps of the algorithm. 
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1.7 





(i) 



(ii) 



Figure 1.6: Given (i) a triangulated 2D space mesh, our algorithm constructs (ii) an unstructured tetrahedral spacetime 
mesh. Time increases upwards in (ii). In this example, the wavespeed at any point in spacetime is one of two 
distinct values: the maximum wavespeed occurs inside a circular cone where the tentpoles are shortest, the minimum 
wavespeed occurs everywhere else. The size of spacetime elements adapts to both changing wavespeed to maintain 
causality as well as to a simulated error metric which depends on the temporal aspect ratio. 



The original Tent Pitcher algorithm proposed by Ongor and Sheffer JUSOOl) applied to one- and two-dimensional 
space domains. The algorithm could guarantee progress only if the input triangulation contained only angles less 
than 90 degrees and if the wavespeed did not increase or only increased smoothly. Ungor and Sheffer referred to 
the causality constraint as an angle constraint and assumed that this angle constraint was Lipschitz-continuous. They 
implemented Tent Pitcher and investigated several different heuristics to improve its performance. 

Erickson et al. (|EGSU02| | extended Tent Pitcher to arbitrary spatial domains, even those with obtuse angles, in 
arbitrary dimensions by imposing additional constraints, called progress constraints. The progress constraint limits 
the amount of progress in time when some vertex of the simplex is pitched, in addition to the limit imposed by causality. 
The progress constraint is a function of the shape of the simplex and is necessary to guarantee progress only when the 
initial space mesh contains a non-acute angle. 

The Tent Pitcher algorithm due to Ungor and Sheffer dUSOOl |U02| | and extended by Erickson et al. ( iEGSU02l) 
applied to the case where the wavespeed at a given point is either constant, decreasing, or increasing smoothly as 
a Lipschitz function. When the wavespeed changes, the previous algorithms take the global upper bound on the 
wavespeed and use that as a conservative upper bound on the wavespeed at every point. Limiting the progress at each 
step by a function of the global maximum wavespeed is unnecessarily restrictive. One would like an algorithm that 
adapts to increasing wavespeeds so that fewer spacetime elements, and therefore less computation time, are required 
to mesh a given volume. We develop such an algorithm, an extension to Tent Pitcher, in ChapterlS] 

The size of spacetime elements affects the numerical accuracy of the approximate solution. Where the exact 
solution is constant or changing very little, large elements suffice to capture this small variation. Parts of the domain 
where the solution changes a lot need to be meshed with smaller elements. In Chapter HI we extend Tent Pitcher to 
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Figure 1.7: The first few tents pitched by Tent Pitcher (lEGSU02b . 

make the size of spacetime elements adaptive to a posteriori error estimates. 

Tent Pitcher is a good alternative to standard time-stepping methods. Global time-marching schemes advance the 
solution from one time value to another everywhere in the domain at once. Computing the solution for the new time 
value in this way has at least one of two disadvantages — either the maximum possible time step is constrained by the 
worst-quality element in the space mesh, or an arbitrarily large coupled system has to be solved simultaneously. The 
former would correspond to the constraint that the height of each tentpole erected by Tent Pitcher must be the same, 
and the latter would occur if the causality constraint were not satisfied by Tent Pitcher. Thus, the advancing front 
approach avoids both limitations associated with time-marching schemes. (Some recent advances have been made on 
adaptive explicit time-stepping methods for ODEs (IEJL03I) .) 

Chapter summary 

In this chapter, we gave a concise introduction to the properties of hyperbolic partial differential equations (PDEs) as 
they concern this thesis. We briefly discussed the challenging problem of generating good quality meshes to support 
fast and accurate numerical solutions to PDEs. We pointed out the advantages of spacetime discontinuous Galerkin 
(SDG) methods over conventional finite element methods. For more about DG methods, the reader is referred to the 
survey ( ICKSOOI l of discontinuous Galerkin methods and their applications. 

We outlined the specific problem of meshing in spacetime to support SDG methods for solving hyperbolic PDEs, 
such as those modeling wave propagation. Meshing directly in spacetime presents unique challenges not encountered 
in classical meshing problems. It is hard to formulate geometric criteria for a good quality spacetime element and a 
good quality spacetime mesh, a problem that is already complicated and appUcation-dependent in Euclidean domains. 
SDG methods, due to their discontinuous formulation, permit nonconforming meshes and more general mesh adap- 
tivity operations. We will exploit this increased flexibility in developing meshing algorithms, extensions to the Tent 
Pitcher algorithm, in the remaining chapters. 

An advancing front algorithm for generating a spacetime mesh, as exemplified by Tent Pitcher ( |U02[ IEGSU05I i. 
has several benefits. Elements are added to the mesh in small subsets called patches. Each element in a patch has at 
least one outflow facet. Since the front at every step is causal, the front represents all the information that is required 
for future iterations of the algorithm. The elements in a patch are coupled but two different patches erected over the 
same front are independent. Thus, the solution procedure is highly parallelizable. Solving each patch is equivalent to 
performing a small finite element computation limited to the patch. Memory requirements are low because patches are 
created in the order of causal dependence and hence can be discarded when they are no longer adjacent to the current 
front. The computation time required for the simulation is reduced due to the limited dependency only between 
elements in the same small patch and because the number of patches is bounded. 
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Chapter 2 

Basic advancing front meshing algorithm 



Tent Pitcher was the first provably correct spacetime meshing algorithm to support an efficient patch-wise solution 
strategy. Tent Pitcher was proposed by Ungor and Sheffer JUSOOl) and later extended by Erickson et al. (IEGSU02[ 
IEGSU05I i. Our advancing front meshing algorithms are extensions to and augmentations of Tent Pitcher. In this chap- 
ter, we give an alternate derivation and proof of correctness of the Tent Pitcher algorithm. The style of the derivation 
in this chapter anticipates extensions to nonlinear problems and mesh adaptivity, problems that are considered in later 
chapters. Also, the alternate derivation in this chapter fixes an error in the paper by Erickson et al. (|EGSU02| | for the 
case of obtuse front triangles, an oversight that was corrected in a subsequent journal paper dEGSUOSl ). 

Our objective in this chapter is to set up a general framework for advancing front spacetime meshing algorithms and 
highlight Tent Pitcher as one specific member of this class of algorithms. The more sophisticated meshing algorithms 
described in later chapters fit in this general framework. 

The Tent Pitcher algorithm is the first instance of an advancing front algorithm to mesh directly in spacetime. The 
algorithm proceeds by advancing a local neighborhood of the front in every step. The algorithm is simple because the 
geometric constraints that limit the amount of progress at each step are linear or quadratic functions of the coordinates. 
For linear problems, the wavespeed O) is a fixed constant, determined by material properties, or can be conservatively 
bounded by the global maximum wavespeed. When the wavespeed is constant, the amount of progress made by the 
front in each local advancing step is limited only by local constraints. 

2.1 Problem statement 

The input is the initial front To, which is a triangulation of the space domain at time f = 0, and the causal slope a 
(equivalently, the wavespeed co = l/cJ), which is determined by the initial conditions of the PDE. In this chapter, we 
assume that the wavespeed CO is a constant, either because the underlying hyperbolic PDE is linear or because (O is the 
global maximum wavespeed. We give an advancing front algorithm such that for every T G M-" there exists a finite 
integer A: > such that the front Ta- after the kth iteration of the algorithm satisfies T^ > T, which means that the entire 
front has achieved or passed time T. The volume between each consecutive pair of fronts T,- and T,+i is triangulated to 
give a patch. 

We will generate meshes with geometrically non-degenerate elements. Specifically, we prove that each element 
has temporal aspect ratio bounded from below. Recall from Chapter [T| that the temporal aspect ratio of an element is 
the ratio of its height to its duration. 

We say that a front T is valid if for every T G M-" there exists a finite sequence of fronts T, Ti, T2, . . ., T/, where 
Til > T, each front in the sequence obtained from the previous front by advancing some local neighborhood. Therefore, 
our problem can be defined as that of constructing a succession of valid fronts such that the volume between successive 
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fronts is triangulated with a small number of spacetime elements (simplices) each with bounded temporal aspect ratio. 

Definition 2.1 (Valid front). A front x is valid if there exists a positive real Af bounded away from zero such that for 
every T G 5R- there exists a finite sequence of fronts T, T\, T2, . . ., Tk where T^ ^ T, each front in the sequence obtained 
from the previous front by advancing some local neighborhood by At. 

The Tent Pitcher algorithm due to Ungor and Sheffer, and Erickson et al., is obtained by making the following 
assumptions: 

1. The front T = 0, or more generally any constant time function T = T, is valid. 

2. The new front t' G advance(T) is obtained by advancing a single vertex of T by a nonnegative amount, i.e., 
A^ = star(P) for some vertex P of T. 

We say that a front t' is obtained by advancing a vertex P of T by Af > to the vertex P' of the front t' if 
t'(/?) = x{p) + Af, for every other vertex q^ pwt have T'{q) — T{q), and every simplex incident on P' is in one-to-one 
correspondence with a simplex of T incident on P. For any front T, vertex p, and real Af > 0, let t' — advance(T,p, Af) 
denote the front obtained from T by advancing phy At. 

We say that a front t' is obtained by advancing a neighborhood A^ of T by Af > to the neighborhood A^' of the 
front t' if Af = max^^(p^(-p))g^{T'(/?) — t(/?)} and if t\A^ and x'\N' coincide. 

Our solution For one-dimensional space domains, we prove that every causal front is valid. In higher dimen- 
sions, we define progressive fronts, and we prove that if a front is progressive, then it is valid. We give an algorithm, 
a modification of Tent Pitcher, that given any progressive front T, constructs a next front T,+i such that T,+i is pro- 
gressive. The volume between T, and T,+ i is partitioned into simplices. The new front T,+i is obtained by advancing 
a local neighborhood of T, by a positive amount bounded away from zero. The algorithm can be parallelized in a 
straightforward manner to solve several patches simultaneously by lifting any independent set of neighborhoods in 
parallel. Whenever the algorithm chooses to lift a local minimum vertex, it is guaranteed to be able to lift it by at least 
Train > where Tmin is a function of the input and bounded away from zero. 

In this chapter, we will restrict front advancing operations to tent pitching, i.e., advancing A^ = star(/5), for some 
vertex p, to A^' such that ^A^ = dN' = lmk{p). Since pitching a vertex does not change the spatial projection, the 
triangulation of every front is isomorphic to the initial space mesh. In this sense, we say that the meshing algorithm in 
this chapter is nonadaptive because the spatial projection of every front is the initial triangulation of the space domain. 
In this chapter, we will assume that the underlying PDF is linear and that a denotes the causal slope, the reciprocal of 
the global wavespeed, a constant throughout spacetime. Alternatively, if the wavespeed is not constant, let ff denote 
(Tmin, a lower bound on the slope of any cone of influence in any spatial direction. 

2.2 Meshing in ID x Time 

We begin by describing the linear nonadaptive Tent Pitcher algorithm to construct spacetime meshes over one- 
dimensional space domains. The input space mesh, denoted by M, is a one-dimensional simplicial complex, i.e., 
a set of vertices with pairs of vertices, say p and q, connected by a segment pq. Let \pq\ denote the length of the 
segment pq, i.e., the Fuclidean distance between p and q. 



21 



Input: A one-dimensional space mesh M C E' 

Output: A triangular mesh H of M x [0,oo) 

The initial front To is M x {0}, corresponding to time f = everywhere in space. 

Repeat for / = 0, 1 , 2, . . .: 

1. Advance in time an arbitrary local minimum vertex P = {p,Ti{p)) of the current front T,- to P' = 

(p,T,+i(p)) such that Tj+i is causal and Tj^i{p) is maximized. 

2. Partition the spacetime volume between t and T,+ 1 into a patch of triangles, each sharing the tentpole 
edge PP'. 

3. Call the numerical solver to compute the solution everywhere in the spacetime volume between T,- 
and Tj+i as well as the a posteriori error estimate. The solution on T, is the inflow information to the 
solver 

Figure 2.1: Advancing front algorithm in IDxTime. 

At each step /, our algorithm chooses to advance an arbitrary local minimum p of the causal front T, to construct 
the new causal front T,+i such that T;+i(/?) is maximized. Since every front has at least one local minimum vertex, 
e.g., the global minimum, the algorithm has a nonempty subset of candidate vertices to pitch at each step. 

Let AB be an arbitrary segment of the front Tj+i . Then, AB is causal if and only if the gradient of the time function 
T,+i restricted to afiab is less than the slope (j{AB), i.e., if and only if 

F T^i+ilabW — r^, < (y{AB). (2.1) 

The new time coordinate t'(;?) of p is constrained because every segment P'Q of the new front is constrained as in 



Equation 2.1 such that t'(/?) < t(^) + \pq\cT for every segment pq e star(p). 



Note that the inequality in Equation 2.1 is a strict inequality, so the set of all T,+i(;?), such that the front T,+i is 



causal, is an open, convex set, whose supremum value is not feasible. Therefore, to maximize Tj+i(p), we compute 



the supremum time value Tjup that satisfies the causality constraint of Equation 2.1 and then take any time value less 



than Jsup but close enough to T^^p not to cause numerical errors due to loss of precision. 



Figure 2.1 describes our linear nonadaptive advancing front meshing algorithm in IDxTime. The solution ev- 
erywhere in spacetime is computed patch-by-patch in an order consistent with the partial order of dependence until a 
particular target time T >0 has been met. 

In the parallel setting, we repeatedly choose an independent set of local minima of the front T,, equal to the number 
of processors, to be advanced in time simultaneously. The resulting patches can be solved independently. If a patch is 
accepted, the local neighborhood of the front is advanced without any conflicts with other patches. 

It is easy to see that each spacetime triangle has one causal inflow facet (a segment) on the old front T; and one 
causal outflow facet on the new front t,+i. 

2.2.1 Progress guarantee 

It remains to show that the algorithm creates only a finite number of triangles to mesh the volume M x [0, T] for any 
target time T > and that the minimum temporal aspect ratio of each triangle is bounded. Both requirements follow 
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from the following lower bound on the height of each tentpole. Recall from Section 1.2.1 that, for a vertex P of 
the one-dimensional front T, the distance Wp is the minimum distance of p from its nearest neighbor in the spatial 
projection of T. 

Theorem 2.2. Let X be a causal front and let p be an arbitrary local minimum of T. Then, for every Af such that 
Q < ht < WpO, the front t' = advance{T,p,At) is causal. 



Proof of Theorem \2. 2\ Only the segments of the front incident on P = {p,T{p)) advance along with p. Consider an 
arbitrary segment pq incident on p. Since p is a local minimum, we have T{q) > T{p). We have 

T\p)<T{p)+At 

< z{p) + Wpa 

<T{q) + \pq\aiP'Q) 



Therefore, the slope of the segment P'Q is less than a{P'Q) and hence P'Q is causal. 



n 



Since each front T is causal, the difference between the maximum and the minimum time coordinate of points of T 
must be smaller than the diameter of the space domain M times the slope a due to causality. Therefore, the algorithm 
must eventually advance the global minimum vertex of T. Thus, we have the following theorem. 

Theorem 2.3. For every i > 0, if the front T; is causal then T; is valid. 



Proof of Theorem \2. 3\ Let Wmin denote the smallest length of the spatial projection of a segment on any front. Since 
causahty limits the gradient of the front, every vertex is pitched eventually when it becomes a local minimum; specif- 
ically, a vertex becomes a global minimum after only a finite number of iterations of the advancing front algorithm. 



By Theorem 2.2 each iteration advances the front in time by a bounded finite amount. Therefore, for any target time 

r > 0, after a finite number of iterations, the entire front has advanced past T . Specifically, we show next that after at 

most 

n(r + diam(M)(Jniax) 

WmmO 

iterations, the entire front has advanced up to or beyond the target time T . Here, n denotes the maximum number of 
vertices and Wmin denotes the minimum spatial length of any segment on any front constructed by our algorithm. Let 

?mm = l^min<7. 

Consider step /+ 1 of the algorithm. By Theorem 2.2 the front Tj+i such that Ti{p) < T,+i(/:i) < Ti{p) + Tj^i^ is 
causal. Therefore, we have shown that if T, is causal then there is a front T,+i = advance (T,,/:i,Af) such that T,+i is 
causal for every At e [0, Tmin)- Note that LpGV(M) '^i+i{p) — ^min + llpev{M) '^i{p)- By induction on /, and because a is 
finite and M is bounded, there exists a finite k > i such that the front ^ satisfies 

Y, T:k{p) > n(r + diam(M)a) 

pGV(M) 



for any real T. Hence, 



max Tii{p) > - 
pev{M) n 



\peviM) , 



>r + diam(M)c7 
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Since Tj^ is causal 



I max Tkip) I — I min Tifp) 1 < diam(M)<7. 

\pSV{M) ) \pev{M) ' 



Since maXp^y/j^x Tii{p) > T + diam{M)a, it follows that minpg^/(^^ T^(/:') > T and so T; is valid. D 

In general, the Tent Pitcher algorithm is free to pitch any vertex, not just local minima; however, the progress 
guarantee of Theorem IZ2] applies only to local minima. Any causal front can be input as the initial front Tq; thus, the 
solution can be saved after each step and the algorithm can be restarted from the last front. 

2.3 Meshing in 2D x Time 

In this section, we consider our advancing front algorithm in 2D x Time. 

2.3.1 Need for progress constraints 

It has been observed by Ungor and Sheffer and also by Erickson et al. that in spatial dimension d >2 the causality 
constraint is not enough to guarantee that spacetime elements created by Tent Pitcher are non-degenerate. With 
causality constraints alone, it was shown by Ungor and Sheffer jUSOOl l, and by Erickson et al. JEGSUOZl l that if 
the space mesh contains an obtuse or a right triangle then Tent Pitcher will eventually construct a front such that no 
further progress is possible while maintaining causaUty. This is the case even for linear PDEs where the wavespeed 
everywhere in spacetime is a fixed constant. 

Our algorithms advance a neighborhood A^ of the front T to the neighborhood A^' of a new front t' where dN = dN' . 
It is necessary that the lower-dimensional simplices of dN satisfy gradient constraints stricter that the causal cone 
constraint. 

Erickson et al. introduced so-called pro^re^i constraints on the front at each stage. Causality limits the magnitude 
of the gradient of every simplex on each front. The progress constraint is a gradient constraint on certain lower 
dimensional faces. 

The progress constraint imposed on T is a gradient constraint on certain edges of the front T. For every local 



minimum vertex p of the front T, the progress constraint limits the gradient of the edges of link(/?). Figure 2.2 is the 
vector diagram that indicates the progress constraint imposed on a single triangle of the front T. 



Interpreting the vector diagram Consider a single triangle PQR of the front T. Let pqr be the spatial 

\pqr 



projection of APQR. Triangle PQR defines a plane in spacetime which is the graph of the linear function t| . : E^ 



^2 

ascent along the plane of APQR and its magnitude is the slope of this plane 



. The gradient of t| , denoted by V t| is a vector in M . The direction of V t| is the direction of steepest 



Figure 2.2 shows the spatial projection Apqr on the left and the vector space M on the right. A vector in M can 



be indicated by an arrow with its tail at the origin. In Figure 2.2 we simplify the depiction of vectors and indicate a 



vector by a point representing where the head of the arrow would be. 

Triangle PQR is causal if and only if the point representing t| lies strictly inside the disk centered at the origin 



with radius equal to a{PQR). In Figure 2.2 we indicate the global slope a by a disk centered at the origin with 
radius a. 



24 




Figure 2.2: Vector diagram: Triangle PQR is causal if its gradient vector is inside the circular disc of radius a. To 
ensure progress in the next step, the gradient vector must also lie outside the shaded red region. [Damrong Guoy, 
personal communication] 
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Suppose we advance P to P' giving the triangle P'QR of a new front t'. Lifting P keeps the gradient along affqr 
unchanged, i.e., t'I^^^ = TJ^g , while increasing the magnitude of the component orthogonal to qr in the direction 
of p. 



Vectors OP, OQ, and OR in Figure 2.2 are orthogonal to qr, pr, and pq respectively and oriented in the direction 



of the third vertex in each case. These three normal vectors subdivide M into three zones. Each zone is identified 
with a vertex of Apqr because whenever V t| lies in that zone the corresponding vertex is a local minimum vertex 
of APQR. For instance, if the gradient vector V t| . lies in zone p, bounded by the vectors OQ and OR, then 

"^{p) < t(^) and x{p) < T(r). 



Consider the situation depicted in Figure 2.2 i). Points and 1 indicate respectively the gradient vector of APQR 
where p is a local minimum vertex and the gradient vector of AP'QR after pitching P to P' so that q is the new local 
minimum vertex. When P is pitched to P', the gradient vector advances from to 1 in the direction OP. Advancing 
from point to point 1 represents positive progress because both points are inside the disk (the corresponding triangles 
are causal) and p is no longer a local minimum. 



Figure 2.2 ii) indicates the possible situation in the next iteration of the algorithm when the new local minimum Q 



is advanced — the gradient vector makes positive progress along direction OQ from point 1 to point 2. If the algorithm 



is too greedy and advances from point to point 3 in the current step, as shown in Figure 2.2 iii), then sufficient 
progress cannot be guaranteed in the next step because the gradient vector may leave the disk (point 4) and violate 
causality. Therefore, it is necessary to ensure the gradient vector of AP'QR does not enter the shaded forbidden zone 



shown in Figure 2.2 iii). The forbidden zone in Figure [Z2[ iii) limits the magnitude of the gradient of the edge pr on 
the new front t', i.e., it imposes an upper bound on T'(p) — T(r). 

Symmetrically, when the local minimum q is pitched from Q to Q' , there is a corresponding forbidden zone that 



limits the slope of the edge Q'R. See Figure 2.2 iv). 

Note that the forbidden zones (and therefore the progress constraint) depends solely on the shape and orientation 
of Apqr and is invariant with respect to uniform scaling. 

2.3.2 Progressive fronts in 2D x Time 

Advancing a vertex p in time increases the gradient of simplices in star(/:') but leaves unchanged the gradient of the 
lower-dimensional simplices in \mk{p). For each face / e link(/:i), the gradient of the front in a directional orthogonal 
to / increases. Our algorithm advances an arbitrary local minimum vertex p of each front T to get a new front t'. The 
new front t' after advancing p is causal if and only if the gradient of / is strictly less than the slope of the causal cone 
constraint that limits the gradient of the simplex F D f such that F e stai{p). 

Lemma 2.4. Let APQR be a triangle of the causal front T with P as its lowest vertex. Without loss of generality, 
assume t{p) < x{q) < T(r). Let (7 > be a slope such that O < l/a){PQR) and let dp denote dist(;?,aff^r). Let £ be 
any real number in the range < £ < 1. Consider two cases: 

L Case /.pqr < n/2: //'||V t| || < (1 — e)(J, then the front x' = advance{T,p,At) satisfies ||Vt'|| < a for every 
At € [0,edpa]. 

2. Case %I2 < ^pqr < tt; /f ||V t| || < (1 —£)a sin Zpqr, then the front t' ~ advance {t , p , At) satisfies ||Vt'|| < a 
for every At £ [0,edpa]. 
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(a) (b) (c) 

Figure 2.3: Triangle pqr of the front T where x{p) < x{q) < T(r) 



Proof of Lemmau^ Let p be the orthogonal projection of p onto the line aff ^r. 

As long as T'(p) < T{q), we have ||Vt'|| < ||Vt|| and the statement of the lemma follows. Hence, for the rest of 
this proof, assume T'{p) > i{q). 

Let Uqr denote the unit vector normal to qr such that n^r ■ (p — q) > 0. Let Vcy,- be the unit vector parallel to qr such 
that ygr • (r — q) > 0. Then, {n^r, v^,} form a basis for the vector space M^. 

Lifting p does not change the gradient of the time function restricted to the opposite edge, so V t' • v^,, = V t • v^y, , 
i.e., V t'I — V t| . The scalar product V t' • n^r can be written as 



Vr'-n^, = 



t'(p)-t(/j) 



\PP\ 



Since q is the lowest vertex of qr, we have Vt' • v^, > 0. Also, we are given that Vt- v^, < (1 — e)c7 < a. Note that 
||Vt|| < Vr-Vf/r- Therefore, ||Vt'|| < a if and only if 





\pp\ 


iiv.gp 



(2.2) 



Case 1: Zpqr is non-obtuse. See Figure 2.3 b)-(c). In this case, we have 'z{p) > T{q) > t{p). We have 

'^'(p) = T(/:') + Ar <T(p)+Ar <T(p) + ea|/7^| 
Since < e < 1 we have e < Vl — (1 — e)^. Therefore, 

t'{p) < X{p) + Vl - (1 - e)2 G\pp\ = T(p) + VCJ^ - (1 - 6)2^2 |pp| 

Since ||V t|^^|1 < (1 -e)(T, we have 



t'(p)<t(p) + |pp|Vct2-|1Vt|^J| 



which is precisely the causality constraint of Equation 2.2 



Case 2: Zpqr is obtuse. See Figure 2.3 a). In this case, we have x{p) < x{q). 
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Let j3 = |p<?|/|pp|. Since \pq\ ^ 0, we have 

X'{p)-X{p) ^ t'{p)-x{q) ^ T{q)-T{p) \pq\ ^ x' [p) - x{q) 

\pp\ \pp\ \pq\ \pp\ \pp\ 

Using Equation |2. 3 1 the causality constraint (Equation|2.2|l can be rewritten as 



-j3||VT 



qrw 



(2.3) 



^M_iM < y^^q^^_^i,VTij 



(2.4) 



We are given that \\V tL^|| < (1 — e)GsinZpqr. Substituting this upper bound on ||V t| || into Equation 



2.4 



obtain the following constraint: 

^ (P) - <g) ^ ^ fVl-{l-efsm^Zpqr) - cj(l - e)j3 sinZp^r 



\PP\ 



(2.5) 



which impUes the causality constraint of Equation 2.4 Since T{p) < x{q) and T'{p) < x{p) + eG\pp\, we have 



t'(p) - T:{q) ^ t:'{p) - t:{p) ^ e(y\pp\ ^ ^^ 



\PP\ 



\PP\ 



\PP\ 



Equation 2.5 is satisfied if 



or equivalently 



e < V 1 - (1 -e)^sin^Z/9^r- (1 -e))3sinZp^r 



(e + (l -e)j3sinZ/?^r)^ + (l -e)'^sin^Z/7^r < 1 



Let A denote the left-hand side of the last inequahty above. 
We have 



A = (e + (l-£)j3sinZpg'r)^ + (l-e)^sin^Z/7^r 
= e^ + (1 - e)^j3^sin^ Zp^r + 2e(l - e))3 sinZp^r + (1 - e)^ sin^ Apqr 
= e^ + (1 - e)2(j32 + 1) sin2 /Lpqr + 2e{l - e)^ sinZpqr 



Note that 



j32 + l = 



\pq\' + \pp\' IwP 



12 



\ppr \pp\ 

by Pythagoras' theorem applied to Apuq. Also, sin Zpqr — \pp\/\pq\- Hence, 



(j32 + l)sin2Z/?^r=L 



Additionally, 



j3 sin Zpqr - 



\pq\ \pp\ 
\pp\ \pq\ 



cos /Lpqr 
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(a) 



(b) 



(c) 



Figure 2.4: Triangle pqr with (a) /^pqr obtuse, so 0^^ = sinZ/jgr < 1, (j){pqr) = ^q/, (b) /^pqr and ^prq acute, so 
^qr= 1, ^{pqr) — sin Zrpq < 1; and (c) Zprq obtuse, so (j)gr — sin Zqrp < 1, (l){pqr) = 0^,- 



Therefore, 



A = e^ + (l — e)^-2e(l -e)cos Zpqr 
= e^ + 1 + e^ - 2e - 2e(l - e) cos Zpqr 
= l+2e^-2e-2e{l-e)cosZpqr 
= 1 — 2e(— e + 1 + (1 -e) cos Zpqr) 
= I -2e{l-e){l+ cos Zpqr) 
<1 



The last inequality follows because Apqr is non-degenerate, so cos Zpqr ^ — 1; also, both e and (1 — e) are positive. 
Therefore, the causality constraint is satisfied. D 



The progress constraint in 2D x Time is motivated by the Lemma 2.4 which states that if APQR is causal and satis- 
fies progress constraint a = l/(o{PQR) then AP'QR obtained by advancing a lowest vertex P in time by e dist(p, aff ^r) a 
is also causal. 

For each simplex (an edge or a triangle) / in star(/?) Ulink(/?), causality and the progress constraint limit the 
gradient of /, separately for each simplex /. Therefore, without loss of generahty, we consider each such simplex / 
separately while defining and computing gradient constraints on /. 

The progress constraint is parameterized by £ and a wavespeed a. We are free to choose any value for e in the 
allowed range. The bound on the gradient of T due to the progress constraint is proportional to the local geometry of 
the spatial projection of T. 

Definition 2.5. For a triangle pqr, for any arbitrary edge of the triangle, say the edge qr, define (^q,- := 1 if both Zpqr 
and Zqrp are non-obtuse; otherwise let (j)qr denote the sine of the obtuse angle of Apqr, i.e., <j)qr := max {sin Zpqr, 
sin Zqrp } Define <j) (pqr) to be the minimum <pe over every edge e of the triangle. Note that < <j) (pqr) < 1 and 



^(pqr) < 1 if and only if Apqr is obtuse (Figure 2.4). 



Definition 2.6 (Progress constraint a). Fix e in the range < e < 1. Let PQR be an arbitrary triangle of a front T. 
We say that the triangle PQR satisfies progress constraint (7 if and only if for every lowest vertex P, equivalently for 
every highest edge QR, we have 



IV tI 



qr\ 



T:{r) - T:{q) 

\qr\ 



<(l-e)(j(^, 



qr- 
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Input: A bounded interval M cM. 

Fix e e (0, 1), a parameter to the algorithm 

Output: A triangular mesh ii of M x [0,oo) 

The initial front To is M x {0}, corresponding to time f = everywhere in space. 

Repeat the following sequence of steps for / = 0, 1 , 2, . . .: 

1. Advance in time an arbitrary local minimum vertex P = {p,'^i{p)) of the current front T, to P' — 
{p, T,+i (p)) such that T,+i is progressive and T,+i (p) is maximized. 

2. Partition the spacetime volume between T and T,+i into a patch of triangles, each sharing the tentpole 
edgePP'. 

3. Call the numerical solver to compute the solution everywhere in the spacetime volume between T, 
and T,+i as well as the a posteriori error estimate. The solution on T, is the inflow information to the 
solver and the solution on T,+i is the outflow information. 

Figure 2.5: Advancing front algorithm in 2D x Time. 



We say that a triangle PQR is progressive if it satisfies progress constraint a. We say that a front T is progressive 
if every triangle of T is progressive. Note that every progressive triangle or front is also causal. 

Suppose a lowest vertex P is being advanced. As long as P is the lowest vertex of APQR, the progress con- 
straint limits II V t'I II but ||V t'L,.|| = ||V t| ||. Assume without loss of generality that T{p) < T{q) < T(r). When 
t'(/?) > t(^), the new lowest vertex is q, so the progress constraint Umits ||V t'| ||. (We can interpret the progress 
constraint inductively as a causality constraint on the 1 -dimensional facet pr opposite q where the relevant slope is 
(l-e)c70,p.) 

In the parallel setting, we repeatedly choose an independent set of local minima of the front T,-, equal to the number 
of processors, to be advanced in time simultaneously. The resulting patches can be solved independently. If a patch is 
accepted, the local neighborhood of the front is advanced without any conflicts with other patches. 

In subsequent sections, we will successively strengthen the progress constraint and correspondingly strengthen our 
notion of progressive fronts. When we refer to a front T simply as a progressive front, we will implicitly understand 
this to mean that T is causal and T satisfies the progress constraint as defined in the most recent definition. 

The value of Ti^i{p) is constrained from above by causality and progress constraints separately for each of the 
simplices incident on p. The final value chosen by the algorithm must satisfy the constraints for each such triangle. 
Therefore, it is sufficient to consider each triangle pqr incident on p separately while deriving the causaUty and 
progress constraints that apply while pitching p. 



The progress of the front at the (/+l)thstep is defined to be T,+i(/7) — T,(/:i). We have already shown in Lemma 2.4 
a lower bound on the amount of progress when subject to causality constraints alone. Next, we prove a lower bound on 
the amount of progress when subject to progress constraints alone. Thus, we obtain a lower bound on the worst-case 
amount of progress when subject to all causality and progress constraints simultaneously. 

The next lemma states that if APQR is causal and satisfies progress constraint a = l/co{PQR) then AP'QR 
obtained by advancing a lowest vertex P in time by (1 — e) ^{pqr) dist(/:i, aff qr) a also satisfies progress constraint a. 
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Lemma 2.7. Let PQR be a triangle of the causal front T with P as its lowest vertex. Let (7 > denote a slope and 
let dp denote dist(/?, aff qr). Let £ be any real number in the range < £ < 1. 

If APQR satisfies progress constraint (7, then for every At G [0, (1 — £)(7dp] the front t' ~ advance (t,pq, At) 
satisfies progress constraint (7. 



Proof of Lemma^y\ By Definition 2.6 the front t' satisfies progress constraint a if and only if 



Vt'I,J| < (l-e)</)p,(T 



!;»/ 



Vt'I,I| < (l-e)V(^ 



\cir 



Vt'I.J| < (l-e)0„a 



\rp 



Since APQR satisfies progress constraint a, we have 

||Vt|,,,J| < {l-e^p,a 
\\Vt\J < (l-e)0„-a 

||VT|,p|| < {l-£)(j,rpCJ 

Also, p is a lowest vertex of t; so T{p) < min{T{q),T{r)}. 

Since advancing p in time does not change the time function along QT, we have ||V t'| || = ||V t| || < (1 — e)0^r<7- 
Consider the constraint || V t'I^^H < (1 - £)^pr(y- As long as T'(p) < t(^), we have || V t'I^^H < || V t|^j^|| < 

(1 — £)^pqG. Similarly, as long as T'{p) < T(r), the constraint || V t'|, J| < (1 — £)(prp(y is automatically satisfied 

because ||VT'|,p|| < ||VT|,p|| < (l-e)0,;,(T. 

When t'(;?) > t(^), the constraint || V t'I^^H < (1 - £)<^pq(y is equivalent to T'{p) < x{q) + \pq\{l - e)(j)pc,a. We 

have 

T'{p)<T:{p) + {l-e)adp 

< x{q) + (1 - £)amm{\pq\,\pr\}<^qr 
<x{q) + {\-£)^pqO\pq\ 

<T{q) + {\-£)^{pqr)a\pq\ 

The last inequality follows because (j){pqr) < (j)pg. 

Similai'ly, whenT'(/:i) > T(r), the constraint ||V t'| || < (1 —e)^^^^ is equivalent to T'{p) < T{r) + {l — £)^rp(y\pr\- 
We have 

T:'{p)<T:{p) + {l-e)(jdp 

< T{r) + (1 - e)amm{\pq\,\pr\}(l)rp 
<T(r) + (l-e)0,pC7|r/?| 

<T{r) + {l-e)^{pqr)a\rp\ 

a 



Note that we have in fact proved the following two statements that imply Lemmas 2.4 and 2.7 
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Let APQR be a triangle of the front T with P as one of its lowest vertices. Without loss of generality, assume 
'^{p) < t(^) < T(r). Let e be any real number in the range < e < L Let a > be any slope. 

1. If ||V t| II < (1 - e)a(j)gr, then AP'QR satisfies || Vt'|| < a for every x'{p) in the range T{q) < x'{p) < x{q) + 



\qr\ 



e<7dist(/7,affg'r). 
2. AP'QTJ satisfies 



IVt'I l|<(l-e)cT^ 



■pq 



and 



|Vt'| ||<(l-e)c7^ 



'pr 



for every 'l'{p) in the range z{q) < T'{p) < T{q) + (1 — e)(Jmin{|/7q'|, \pr\]. 

The above claims are stronger than Lemmas |2.4| and [ZTJ because we do not assume that APQR and the front T is 
progressive. 



As a corollary to Lemmas 2.4 and 2.7 we conclude a positive lower bound on the minimum tentpole height, 



similar to that obtained by Erickson et al. ( |EGSU02| l in the nonadaptive case when the wavespeed everywhere is a 
fixed constant co. Let a denote the causal slope. 

Corollary 2.8. Let X be a causal front. Let e be any real number in the range < £ < 1. If every triangle APQR 
of 1 satisfies progress constraint O, then, for any local minimum p of 1, the front x' — advance(x,p^ht) is causal and 
satisfies the progress constraint O for every Af £ [0,min{e, (1 — e)}wp(7]. 

Proof of Corollary \2. 8\ Causality and progress constraints apply to each triangle PQR of the front incident on P. 
We show that each constraint individually guarantees a positive lower bound of min{£,(l — £)}wpff on the height 
i'{p) ~ t{p) of the height of the tentpole at p. Therefore, all the constraints together ensure a tentpole height of at 
least as much. 

to each triangle PQR, we conclude that the front t' is causal for every At e [0, ewpO]. 



Applying Lemma 



2.4 



Consider an arbitrary triangle PQR incident on P. Advancing p in time does not change the gradient of the front 
along qr. Therefore, the gradient of the new front is limited only along pq and pr by the progress constraint. By 
Lemma [24} we conclude that each of the gradient constraints along pq and pr guarantee a tentpole height of at least 



(1 — e) dist(/:i, aff^r) a. Applying Lemma 2.4 to each triangle PQR, we conclude that the front t' satisfies the progress 
constraint a for every Af e [0, (1 — e)wpa]. 



U 



Theorem 2.9. For any i > 0, if the front Tj is progressive then T,- is valid. 



Proof of Theorem 2.9 The proof is identical to that of Theorem 2.3 and follows because we choose £ > so that 
Tmin > is bounded away from zero. D 

We thus have the following theorem. 

Theorem 2.10. Given a triangulation M of a bounded planar space domain where Wmin is the minimum width of a 
simplex of M and <J is the minimum slope anywhere in M x [0, °°), for every £ such that < £ < 1 our algorithm 



constructs a simplicial mesh of M x [0, T\ consisting of at most 



n(T +diam(M)a) 
min{£,l — £}o"w,„/, 



spacetime elements for every real 



T >0 in constant computation time per element, where A is the maximum vertex degree. 

Proof of Theorem \2. 1 0\ By Lemmas |2.4| and [ZT] it follo ws t hat the height of each tentpole constructed by the algo- 
rithm is at least T^^ = min{£, 1 — £}ffWmin- By Theorem 



2.9 



after constructing at most k < 



n(7'+diam(A?)g) 
T 



patches. 



the entire front Xj^ is past the target time T. Since each patch consists of at most A elements, where A is the maximum 
number of simplices in the star of any vertex of M, the theorem follows. D 
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2.4 Meshing in arbitrary dimensions dD x Time 



For spatial dimension li > 3, we prove results analogous to those in Section 2.3 For each front T,, causality limits the 
gradient of 1 -dimensional simplices, and causality and progress constraints together limit the gradient of ^-dimensional 
simplices for each 2 <k <d. 

Pitching vertex po, a local minimum of the front T, alters the gradient of all faces of star(/7o) of all dimensions 
and each of these faces limits the progress of po in time because of causality and progress constraints on their gradient 
in time. Our objective is to impose gradient constraints on each ^-dimensional simplex such that each constraint 
individually guarantees a minimum positive progress in time for po- Hence, all these gradient constraints taken together 
also imply positive progress which in turn implies a lower bound on the height and therefore the temporal aspect ratio 
of each spacetime element. 

Fix e G (0, 1). The progress constraint is parameterized by e and the wavespeed a. We are free to choose any 
value for £ in the allowed range. 

Consider an arbitrary A:-dimensional face foA^2- ■ -Pk incident on a vertex Pq of a front T where 2 < k < d. Let 
PoPiPi- ■ -Pk denote the spatial projection of this simplex. The extent to which progress constraint a limits the gradient 
ofPoPiP2. ■ -Pk, ie., the gradient of the time function T restricted to Pq A ^2- -Pk, depends on the geometry of the simplex 
PoPiPi- ■ -Pk in space. 

For each integer / such that < / < A:, let A, denote the facet popiP2- ■ -Pi-iPi+i ■ ■ -Pk opposite /?,, let pi denote the 
orthogonal projection of p, onto aff A,, let m; denote the point of A, closest to /5,, and let 0, denote Zp/M,/?,-. We say that 
the simplex popip2- ■ -Pk is obtuse if /?,■ ^ A,- for some /. 

Definition 2.11. For the k-dimensional simplex poPiPi- ■ -Pk define 

(l){P0PiP2---Pk)-^ min{sin0,} 

0</<A; 



As in 2D x Time, there are two key lemmas in higher dimensions that imply a lower bound on the worst-case height 
of any tentpole. Lemma 2.12 was proved implicitly in the paper by Erickson et al. dEGSUOit but was applied only to 



the case of constant wavespeed. 

The first lemma states that if PoPiP2- ■ Pk is causal and satisfies progress constraint a then PqP\P2. . .Pk obtained by 
advancing a lowest vertex Pq in time by £(jdist(po,aff/?ip2- • -Pk) is also causal. 

Lemma 2.12. Let PqP\P2. . .Pk be any k-dimensional simplex of the causal front T for arbitrary I < k < d with Pq as 
its lowest vertex. Let po denote the orthogonal projection of po onto affpip2. . .pk. Let (J > be a slope such that 
G < l/co{PoPiP2...Pk) <^nd let dpg denote dist{po,SLffpip2...Pk) — \poPo\. Let £ be any real number in the range 
< £ < 1. Consider two cases: 

1. Casepo e PiP2---Pk -^11^ '^Ln, n, II < (1 ^e)^. then, for every At G [0,£(ipQ(j], the front t' = advance{T , po, At) 

satisfies \\Vt'\^^^^^^^^ J < a. 

2. Casepo ^PiP2---Pk : Let u denote the point of the facet pip2. . .pk closest to po. //||V t| || < (1 — £)Gsm/lpoupo, 
then, far every At £ [O^edp^o], the front t' — advance (T,po, At) satisfies ||V t'| || < (7. 
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Figure 2.6: An obtuse tetrahedron with po ^ pipi- ■ -Pk- 



Proof of Lemma 2.12 The constraint on PqP\P2- ■ -Pk is equivalent to the following: 

T^(j.o)-T(p-o) y 2_||v^| ^ 

IPOPoI ^^ ll^%P2...pJI 

Without loss of generality, assume t(/7o) < "^{pi) ^ '^(/'2) < ■ ■ ■ < '^{Pk)- 
Case 1: po e PiPi- • -Pk- In this case, we have t{po) > t:{pi) > t:{po)- Hence, 

'^'iPo)-<Po) ^ t:'{po)-t:{po) 



\PQPo\ 



\POPO\ 



Since i'{pq) — '^{po) < £g\pqpo\, it suffices to show that 



Since 11 V t 



I.e., 



„ II < (1 ^ £)f 1 it suffices to show that 

e2<l-(l-e)2 

2e^ < 2e 



which is true whenever < £ < 1 . 
Case 2: po^PiPi Pk- 

Note that u ^ po, hence, |/'o/'o|/|/'om| = sinZpoupo < 1- We consider two sub-cases separately: (a) t{pq) > t(/7o) 



and (b) t(/?o) < '^(po)- See Figure 2.6 
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Note that we could have merged case 2(a) with case 1; however, then the antecedent of the theorem would have 
depended on the current front T. This would have, in turn, complicated the implementation of the algorithm and the 
proof of the subsequent theorems. 
Case 2(a): T(po) > T(po). In this case, since z{po) > ^(po), we have 

t:'{po)-'^{po) ^ t:'{po)-t:{po) 
\poPo\ ~ \popo\ 

Since t'(/7o) — '^(po) < £(y\pQPo\, it suffices to show that 



Since II Vt| || < (1 —e)asinZ/7oMpo < (1 ^e)^. it suffices to show that 

e2<l-(l-e)2 

i.e., 

2e^ < 2e 

which is true whenever < e < 1 . 

Case 2(b): T(po) < T(po). In fact, we will handle this case with a potentially weaker antecedent than the one in the 
statement of the theorem. Let v denote thepoint of aff/?ip2 ■■•/?<: closest to po among all points such that t(v) > t(/5o). 
Since u e pipi- ■ -Pk we have T{p) > t(/?o); hence, \pov\ < \pou\ and \pov\ < \pou\. Hence, 

smZpovpo = n — r ^ n r = smZ/joM/Jo (2.6) 

\pm |pow| 

We will prove the following statement which, by the inequality of Equation |2.6| implies the theorem: 

If 11^ '^lpiP2 PtII < (l^£)o'sinZ;?ov^o, then, for every Are [Ojel/^o^olc], the front t' = advance (T,/?o,Af) 

satisfies ||VT'|p„^,^^...pJ|< (J. 

Let j3 = |/'ov|/|po/'o|- By our choice of v, V t| is a vector pointing in the direction from po to v. Hence, 

'z^'(Po) - </5o) t'(/?o)-t(v) 



bopol IpopoI 



^ll^^lp.P2...pJl (2.7) 



Using Equation 2.7 the constraint on PqPiP2- ■ -Pk can be rewritten as 



Now, 



x'{pQ)-x{y) ^ x'{po)~x{po) 



\PQPO\ \PQPO\ 

Since ^'{po) — T(po) < £<y\poPo\^ it suffices to show that 
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By the antecedent, we have || V t| „ |1 < (1 — £)GsmZpovpo. Hence, it suffices to show that 

£ < V 1 - (1 -e)^sin^Z/?ovpo-j3 (1 -e)smZpovpQ 



I.e., 



I.e., 



I.e., 



Note that 



(e + j3(l-e)sinZ/:>ovpo)^< 1 - (1 -e)^sin-Z/7ov/5o 

e^ + (j3(l -e)sinZpovpo)^ + 2j3£(l -e)sinZ/7ovpo + (l -e)^sin^Z/?ov/?o < 1 

e^ + 2e (1 - e)j3 sinZpovpo + (1 - e)^(j3^ + 1 ) sin^ Z/jqv/jo < 1 



,|2_Ll„„.^„|2 

jS' + l 



bopoP 



IPOPO 



\pov[ 

2 



where the last equality follows by Pythagoras' theorem applied to ApopQV. (See Figure 2.6 ) Also, sinZpoi^'Po = 
\PoPo\/\pov\. Hence, (j3^ + l)sin^Z75ov^o = 1- 
Additionally, 

fl • / - \PQv\ \POPo\ 

j3 sin Zpovpo - 



Ipopol \pov\ 
= cos ZpQvpo 

Therefore, it suffices to show that 

e^ + 2e(l -£) cos Zpovpo + (l-e)^ < 1 
Observe that e^ + 2£ (1 - e) + (1 - e)^ = (e + (1 - e))2 = i. Therefore, 

£^ + 2e(l -e)cosZpovpo + (l -e)^ = 1 -2£(1 -e) (1 - cos Zpovpo) 

Hence, it suffices to show that 

1 -2e(l -e)(l -cos Zpovpo) < 1 

Recall that Zpovpo is acute, so < cos Zpovpo < 1, i.e., 1 — cos Zpovpo > 0. Hence, the claim follows whenever 
0<e<l. n 

2.4.1 Progress constraint for a fc-dimensional simplex 



Motivated by Lemma 2.12 we derive progress constraints that are sufficient to ensure that the antecedent of Lemma 2.12 
is satisfied in the next step by the front t'. 

Consider an arbitrary ^-dimensional face popiP2- • -Pk incident on po for any k in the range 2 < k <d. 
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Definition 2.13 (Progress constraint). We say that a simplex PoPiP2- ■ Pk of the front z satisfies the progress constraint 
(7 if and only if for every highest facet Ai where < i <k we have || V t|^. || < (1 ^ £) <P{poPiP2- ■ -Pk) C- 

Note that if PoPiP2- ■ Pk with Pq as its lowest vertex satisfies the progress constraint a (Definition |2.13| l then 

11^ T|pip,...p, II < (1 - £)<l>ipopiP2- ■ ■Pk)cy. 

Lemma 2.14. Let PqP\P2. . .Pk be any k-dimensional simplex of the causal front T for arbitrary 2 < k < d with Pq as 
its lowest vertex. Let <j) denote 0(/?o/'i/'2- ■ -Pk)- For each i with <i <k, let v; denote f/ie/jo/nf q/aff (A,n Aq) closest 

to PQ. 

IfPQPiP2. ■ -Pk satisfies progress constraint (J, then ft? r every At G [0, (1 — e)0(7|v,7?o|] the front t' — advance {t , pq, At) 
satisfies || V t'|^.|| < (1 — £)(j)(7 for every i such that < i < k. 

Proof of Lemma \2J4\ Note that t'|^ = t|^ ; hence, the statement is trivial for / = 0. 



Since PoA^2- • Pk satisfies progress constraint a, by Definition 2.13 we have || V t|^ || < (1 — e)0C7 for every / 
such that 0<i<k. 

Consider an arbitrary / in the range I < i < k such that A, is a highest facet. 

Note that for every simplex / G link{po), the time function T restricted to aff/ is unchanged by pitching po. 
Therefore, pitching po increases V t|^ in a direction orthogonal to aff A, n Aq while keeping the component of T along 
A, n Ao unchanged. Hence, 

Since ||V t|^.|| > and T'(po) < t{pq) + {I - e)(j)a\vipQ\, we have 



\ViP0\ \ViPQ\ 

Therefore, PqPiP2. . .Pk satisfies the constraint in the statement of the lemma. D 



Note that we have in fact proved the following two statements that imply Lemmas 2.12 and 2.14 

Let PqP\P2. . .Pk be any A:-dimensional simplex with Pq as its lowest vertex. Without loss of generality, assume 

t(po) < '^{p\) ^ '^{pi) 1^ ■ ■ ■ 1^ '^{Pk)- For each / with < i < k, let v,- denote the point of aff(AjnAo) closest to po. 

Let e be any real number in the range < £ < 1 . Let <T > be any slope. 

1- If II V ^\pun-,J ^ (1 - £)<y(t>{PiP2- ■ -Pk), then P^PiPj- • -Pk satisfies || V T'\p^„^„^,„J < O for every x'{po) in 
the range T(pi) < T'(po) < T(pi) + ec7dist(po,affpip2- --w)- ■ 

2. P^PiP2---^/i satisfies 

yi,0<i<k IIVt'I . II <(l-e)(70(A;) 

for every z'{po) in the range T{pi) < z'{po) < t(/?i) + (1 - e)(Tmino<,<A:{|v,/?o|}- 



The above claims are stronger than Lemmas 2.12 and 2.14 because we do not assume that PoP\P2. ■ Pk and the 
front T are progressive. 
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2.4.2 Combining constraints for all dimensions 

Pitching vertex p alters the gradients of all faces of star(/:i) of all dimensions and each of these faces limits the progress 
of p in time. We interpret the progress constraints mentioned in the previous lemmas as gradient constraints on lower 
dimensional faces of the front. These gradient constraints are imposed simultaneously on all facets of dimension 
from 2 through d of the front T, at each step. The previous lemmas guarantee that for each k, where \ <k<d, the 
causality and progress constraint for each ^-dimensional face / ensures positive progress of p so that the corresponding 
face /' on the new front also satisfies the relevant causality and progress constraints. 

Thus, every A:-dimensional highest face A has a gradient constraint imposed on it for every {k+ 1) -dimensional 
highest face Y such that A C F. 

Theorem 2.15. For any i > 0, if the front T,- is progressive then T; is valid. 

Proof of Theorem \2. 1 5\ The proof is identical to that of Theorem 2.3 D 

We thus have the following theorem. 

Theorem 2.16. Given a triangulation M of a bounded planar space domain where Wmin is the minimum width of a 
simplex of M and <J is the minimum slope anywhere in M x [0, °°), for every £ such that < £ < 1 our algorithm 



constructs a simplicial mesh of M x [0,T] consisting of at most 
T >0 in constant computation time per element. 



n(T +diam(M)a) 



spacetime elements for every real 



algorithm is at least T^m = ecwmin- By Theorem 



Proof of Theorem \2. 1 6\ By Lemmas |2. 12| and |2. 14| i t follows that the height of each tentpole constructed by the 

diam(M) C7 



2.15 



after constructing at most k < 



patches, the entire 



front Ti is past the target time T. Since each patch consists of at most A(M) elements, the theorem follows. D 

2.5 Element quality 

Various definitions of quality exist for linear elements in EucUdean space, such as ratio of circumradius to inradius, 
circumradius to shortest edge length, radius of minimum enclosing ball to maximum enclosed ball, etc.. Each defi- 
nition is useful to quantify the suitability of the element for the particular numerical technique being used. In E^, all 
measures of quality are within constant factors of each other, which is not true in higher dimensions. We recommend 
the excellent survey due to Shewchuk (IShe02b on various aspects of element quality. 

The spacetime domain E'' x M is not Euchdean because the time axis can be scaled independently of space. It is 
therefore nontrivial to define the quality of a spacetime element. We consider the temporal aspect ratio of an element 
as a measure of its quality. The temporal aspect ratio of an element was defined in Chapter |2] as the ratio of its height 
to its duration. We prove that our worst-case guarantees on the height of each tentpole lead to a lower bound on the 
worst-case temporal aspect ratio of any spacetime element constructed by our algorithm. 

In this chapter, we have proved a lower bound on the worst-case height of each tentpole, a progress guarantee 
similar to that of Erickson et aJ. (IEGSU02| | and Abedi et al. (lACE+04b . For the linear nonadaptive version of our 
algorithm, where the slope is a everywhere, this progress guarantee can be rephrased as follows: the height of each 
spacetime tetrahedra in the tent pitched at P is at least eowp where e £ (0, j] is a fixed parameter and Wp denotes the 
distance of p from the boundary of the kernel of lmk{p) in the spatial projection. Thus, the height of the tentpole at P 
is bounded from below by a positive function of e, the slope a, and the shape of the triangles in star(/?) . 
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Note that the temporal aspect ratio is always in the range (0, 1] with a larger value corresponding to a "better" 
element. The duration of the tetrahedron P'PQR can be at most 2GWp because both facets PQR and P'QR are causal. 
Together with the lower bound on the height of the tetrahedron, this implies the following theorem. 

Theorem 2.17. The temporal aspect ratio of any tetrahedron in ii is at least e/2. 

When the wavespeed is not constant, the worst-case temporal aspect ratio that we can prove is smaller by a factor 
<7min/<7max, i-C-, proportional to the ratio of maximum to minimum wavespeeds or the Mach factor. 

Bisecting a triangle on the front can improve the temporal aspect ratio of future tetrahedra because the smaller 
triangles may be limited by a smaller wavespeed than their larger ancestors. However, newest vertex bisection does 
not improve the quality of the spatial projection of front triangles, beyond a very limited amount. Repeatedly bisecting 
a single triangle produces smaller triangles from at most four different similarity classes. To improve the temporal 
aspect ratio, mesh smoothing and other generalized mesh adaptivity operations are more useful because they improve 
the quality of the spatial projection of the front. 

2.6 Size optimality 

The computation time of any solution strategy increases with the number of patches — computing the solution within 
each patch is expensive and, in our scenario, much more computationally expensive than the mesh generation, espe- 
cially for high polynomial orders. Therefore, it is important to generate a mesh with close to the fewest number of 
patches necessary for a given accuracy. 

Our algorithm constructs groups of coupled spacetime elements inside each tent such that the boundary of the tent 
is causal. The number of elements in the resulting patch is at most the maximum vertex degree of the underlying 
space mesh. Each element constructed by our algorithm has both a causal inflow facet and a causal outflow facet; 
additionally, the spatial projection of each element PqPqP\P2. . .Pd is the simplex poPiPi- ■ -Pd in the original space 
mesh. Given a triangulation M of the space domain and a target time T, we say that a simplicial spacetime mesh of 
M X [0, r] is solvable if (i) each element has both a causal inflow facet and a causal outflow facet; and (ii) for every 
point X in the spatial projection A of each element, the diameter of A does not exceed the diameter of the simplex of M 
containing x. 

Fix an arbitrary point x in space. The size of a spacetime mesh of M x [0, T] is the maximum over x E M of the 
number of spacetime elements that intersect the temporal segment x x [0, T]. 

Consider the linear nonadaptive instance of our algorithm. The wavespeed everywhere in spacetime is constant, 
equal to a. 

Theorem 2.18. The size of the mesh constructed by our algorithm is (9(l/e ) times the minimum size of any valid 
mesh of the spacetime volume M x [0, T]. 

Proof. Let D and p denote the diameter and inradius respectively of the simplex pop i pi- ■ -Pd of M containing x. 
Causality limits the gradient of any edge of the J-simplex to less than a. Therefore, any temporal segment of duration 
(d—l) oD must intersect at least two distinct elements in a solvable mesh; therefore, the number of spacetime elements 
in a solvable mesh that intersect x x [0, T] is at least YT /{{d — \)(yD)\ . 

Consider a minimal sequence of tent pitching steps, called a superstep, in which each vertex of popipi- ■ -Pd is 
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lifted at least once. When po is pitched, the amount of progress made by x is proportional to dist(x, aff Aq) . Since 

Y^ dist(x, aff A,) > p 

0<i<k 

the amount of progress made by x during a superstep is at least eayp, where 7 e (0, 1] denotes the minimum over 
all / such that < / < A; of vvp./dist(;5,, affA,). Hence, after at most \T /{eGYp)~\ supersteps, the point x achieves or 
exceeds the target time T. 

Causality limits the gradient of any edge of the simplex popipi- ■ -Pd to less than a. Therefore, any d—l vertices 
of pop \P2- ■ -Pd can advance less than 2{d— \)aD units of time total without also advancing the dih vertex. Therefore, 
the number of steps in each superstep is at most \{2{d~ \)aD)/{£aw)\ where w — mino<,<A:Wp.. It follows that 
the number of tetrahedra in the spacetime mesh constructed by our algorithm intersected by jc x [0, T] is at most 
\T/{eayp)}-[{2{d-l)aD)l{eaw)\. 

The ratio of the upper bound to the lower bound on the size is 

^ ^ , 1 D^ 
Old 



e^ ypw 



U 

Note that for every fixed input space mesh M, the size of the spacetime mesh of M x [0, T] constructed by our 
algorithm is at most (9(l/e^) times that of a solvable mesh of the same spacetime volume, where the constant in the 
big-Oh notation depends on M. 

When the wavespeed is not constant, the size-optimality factor that we can prove is larger by a factor (Tmax/Cmin, 
i.e., inversely proportional to the Mach factor. 

When the spacetime error indicator demands different mesh resolution in different portions of the spacetime do- 
main, we can prove that the mesh constructed by our algorithm is still near-optimal where the approximation ratio is 
worse by a factor proportional to the ratio of maximum to minimum element sizes allowed by the error tolerance. 

Our advancing front meshing algorithms are parameterized by real numbers like e that control the extent to which 
they greedily maximize the height (and hence the temporal aspect ratio) of new elements. Being greedy at each step is 
not always a good choice for maximizing the worst-case temporal aspect ratio and for minimizing the total number of 
spacetime elements. We have especially clear evidence that while attempting to satisfy coplanarity constraints in order 
to coarsen parts of the front it may be necessary to reduce the heights of tentpoles. As a result of being less greedy at 



each step, do we actually increase the number of spacetime elements to an unbounded extent? Theorem 2.18 says that 
by being less greedy, our advancing front algorithm is near-optimal in the number of elements compared to a solvable 
mesh. 

2.7 Geometric primitives 

Maximizing the height of each tentpole subject to causality and progress constraints is equivalent to shooting a ray in 
an arrangement of (a small number of) planes and infinite cones. Consider the case where the wavespeed everywhere 
in spacetime is the same, i.e., the slope of every cone of influence is equal to a. 

When a vertex p of the causal front T is advanced in time to get the front t' by pitching a tentpole from P= {p,T{p)) 
to P' = (/7,t'(/7)), the new front t' is causal if and only if for every triangle pqr incident on p the slope of triangle 
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P'QR is less than a. Let Kqr denote the plane through QR making a slope of a with the space domain. The slope of 
AP"QR is equal to a if and only if P" = {p, t"{p)) is on the plane Ttqr and the slope of AP'QR is less than a if and 
only if t{p) < t' {p) < T"{p). Thus, T"{p) can be computed by shooting a ray in spacetime with origin at P oriented 
in the positive time direction and finding the point of intersection of this ray with the plane n^r- The triangle P'QR is 
causal only if \PP'\ < \PP"\. Finding the supremum height of the tentpole at P subject to causality constraint alone is 
equivalent to finding the distance along the tentpole direction before the ray originating at P in the tentpole direction 
first intersects any plane Ttgr among all such planes for every edge qr e link(p). 

The progress constraint ||V t'| || < (1 — e)(j)pqG is satisfied if and only if the edge P'Q is below or on the infinite 
right circular cone Q with apex at 2 = (?, '^{l)) and axis in the positive time direction. Thus, maximizing the height 
of the tentpole PP' is equivalent to shooting a ray with origin at P in the tentpole direction and determining its first 
intersection with Q. 

Thus, finding the supremum height of the tentpole at P subject to causality and progress constraints is equivalent to 
answering ray shooting queries among an arrangement of planes, one for each edge qr G link(/?), and infinite circular 
cones, one for each vertex q neighboring p. 

Running time The running time is 0{deg{p)) for constructing the patch at p, where deg(/?) is the number of 
simplices in star(/5) of any dimension, hence 0(1) amortized per spacetime element. Additional complexity at each 
step comes from the need to choose the local minimum vertex to pitch; the choice of local minimum can be arbitrarily 
complicated. Advancing p destroys at most one local minimum (p itself) and creates at most deg{p) new local minima 
(the neighbors of p). Thus, the set of local minima of each front can be maintained efficiently. 

If the tentpole is inclined, for instance to track a moving boundary or shock interface, then the problem of maxi- 
mizing tentpole height is still equivalent to answering ray shooting queries. The problem is a little more complicated 
if the choice of tentpole direction is not specified and must be chosen to optimize some geometric or numerical quality 
criteria. 

Note that if the wavespeed is not uniform everywhere, the problem is more complicated by the fact that nonlocal 
cones can limit the height of a tentpole and that progress constraints must anticipate future wavespeeds. In problems 
with nonuniform wavespeeds, maximizing the height of a tentpole is not always equivalent to answering a ray shooting 
query — a triangle of the new front can be tangent to a remote cone of influence even if the corresponding tentpole is 
still below the cone. 

We will postpone discussions of geometric primitives for greedily maximizing the height of each tentpole for such 
more general problems to later chapters. 

Chapter summary 

The Tent Pitcher algorithm is the first instance of an advancing front algorithm to mesh directly in spacetime. The 
algorithm is simple because the geometric constraints that limit the height of each tentpole are linear functions of the 
coordinates. The height of each tentpole depends on the size of the space elements incident on the tentpole only and 
not on any global property of the space mesh. Consequently, we were able to prove a lower bound on the temporal 
aspect ratios of spacetime tetrahedra constructed by our algorithm. Progress constraints are artifacts of our algorithm, 
necessitated by its local advancement strategy, and not required explicitly to satisfy causality. Even though progress 
constraints limit the progress in time at each step, we were able to prove that the number of elements produced by our 



41 



algorithm for meshing a given spacetime volume is not much more than that in a size optimal causal mesh of the same 
volume. 

In this chapter, we described the constant wavespeed case where the height of each tentpole is limited only by 
local constraints, as in the case of linear PDEs. We were able to generalize Tent Pitcher to pitch inclined tentpoles 
corresponding to advancing the front in time as well as changing the spatial projection of the front vertices. We were 
also able to incorporate additional operations such as edge flips. These new operations allow the algorithm to recover 
from a bad quaUty initial front (the space mesh M). In Chapter|7j we will also use these operations to track spacetime 
features such as moving boundaries by aligning the tentpole direction and/or causal facets of the front with these 
features. 
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Chapter 3 

Meshing with nonlocal cone constraints 



The basic advancing front algorithm of Chapter |2] relied on a fixed lower bound a ~ Oram on the globally minimum 
slope. The progress constraint was a function of this minimum slope O-iam and this constraint restricted the progress at 
every step of the algorithm. Due to nonlinearity, the speed at which a wave propagates may be different in different 
parts of the domain. Even at the same point in space, the wavespeed may change with time. In general, the wavespeed 
may be a function of the solution. In this chapter, we consider the problem of meshing the spacetime domain when 
causality limits the gradient of mesh facets differently in different parts of the domain because the causal slope is not 
a constant everywhere. 

In this chapter, we make our algorithm respond to changes in the causal slope (the reciprocal of the wavespeed), 
even if such changes are discontinuous and unpredictable, subject to the no-focusing assumption. When the wavespeed 
is not constant, a fast distant wave corresponding to a shallow cone of influence can limit the amount of local progress 
made by our advancing front algorithm. Thus, when the wavespeed is not constant, the most limiting cone constraint 
can be nonlocal. Nonlinear PDEs or even linear PDEs with discontinuous initial conditions can exhibit such behavior 
The no-focusing assumption. Axiom [LT] and its anisotropic version Axiom [T73| allows us to conservatively estimate 
the slope at a point P in spacetime where the solution has not been computed yet, given the cone of influence of every 
point on the current front. The slope o{P) is bounded by the minimum and maximum slopes of all cones of influence 
that contain P. We extend the algorithm of Chapter|2]to construct only causal fronts, making use of the no-focusing 
assumption to estimate the causaUty constraint in the unmeshed and unsolved portion of the spacetime domain ahead 
of the current front. 

3.1 Problem statement 

Just like in Chapter l2] the input to our problem is the initial front Tq and the initial conditions of the underlying 
hyperbolic PDE in the form of the causal slope C7(P) for every point P of Tq. We want an advancing front algorithm 
such that for every nonnegative real T there exists a finite integer k>Q such that the front T^^ after the kih iteration of 
the algorithm satisfies T^ > T . As before, we would like to characterize valid fronts as those that are guaranteed to 
make finite progress at each step. 

The main difficulty in characterizing valid fronts arises when the causal slope at a given point in the space domain 
decreases over time, i.e., the wavespeed increases. For instance, suppose APQR is causal. However, as soon as the 
local minimum vertex, say P, is advanced in time by an arbitrarily small positive amount to P', the triangle P'QR 
may intersect a cone of influence with a much smaller slope, i.e., a{P'QR) <C o{PQR). Consequently, AP'QR is not 
causal. The decrease in causal slope from a{PQR) to a{P'QR) prevents the front from making nontrivial progress by 
advancing the local minimum vertex P. 
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Our solution For one-dimensional space domains, we prove that every causal front is valid. We give an algo- 
rithm that given an arbitrary causal front T; constructs a next front T;+i = advance (T,,/:i,Af) such that T;+i is causal and 
T;-|-i(p) is maximized. 

In higher dimensions, i.e., d > 2, we define progressive fronts and prove that if a front is progressive then it is 
vaUd. We give an algorithm that given any progressive front T, constructs a next front T,+ i — advance (T,,;?,Af) such 
that T,+i is progressive and T,+i(p) is maximized. Whenever p is a local minimum of T,, the progress T,+i(p) — Ti{p) 
is positive and bounded away from zero. 

Our algorithm resolves the following conundrum. The progress of the front at each step ; is limited by the progress 
constraint that must be satisfied by the next front at step / + 1 . However, we do not know what the next front is 
unless we know how much progress is possible at step /. Intuitively, the geometric constraints that apply at any 
given iteration of the algorithm are predicted by simulating the h next tent pitching steps for some parameter h>Q. 
We initially make conservative assumptions about the future wavespeed and successively refine the estimate of the 
wavespeed encountered in the next h iterations. 

Our algorithm is the first algorithm to build spacetime meshes over arbitrary-dimensional triangulated spatial 
domains suitable for solving nonlinear hyperbolic PDEs, where the wavespeed at any point in spacetime depends on 
the solution and cannot be computed in advance. Moreover, the solution can change discontinuously, for instance 
when a shock propagates through the domain. As long as the no-focusing condition (Axiom [TTT] i holds, the resulting 
mesh is efficiently solvable patch-by-patch by SDG methods. 

3.2 Nonlocal cone constraints in ID x Time 

To complete the description of the algorithm in ID x Time, it remains only to describe how to maximize the height 
T;+i (p) — Ti{p) of the tentpole subject to the causality constraint. Even in ID x Time, this optimization problem is not 
trivial in the presence of nonlinearity because of nonlocal cone constraints. Nonlocal cone constraints are handled in 
the same way in higher dimensions but the algorithm is more complicated because of progress constraints. 

Recall that the cone of influence of a point P has its apex at P and its slope in any spatial direction is the maximum 
slope of any characteristic through P in that direction; fast waves correspond to cones with smaller slope. 

When the PDE is nonlinear, a distant but fast wave, i.e., a nonlocal cone, can overtake a slower wave and hence 



limit the duration of new elements. See Figure 3.1 We assume as a conservative estimate that the wave travels through 
the space domain M along the shortest path in the ambient Euclidean space. When the medium is anisotropic, waves 
travel faster in one direction than the other, so the cones are non-circular, for instance, with elliptical cross-sections. 
We will discuss the situation for anisotropic problems in ChapterlTJ it should be noted that the discussion in the current 
chapter is not restricted to circular cones. 

Our algorithms advance a local neighborhood A^ of the front T at each step to the neighborhood A^' of a new causal 
front t'. Hence, t' is causal if and only if (i) the gradient of t' at every point of A^' is less than the minimum slope 
anywhere in A^, and (ii) the neighborhood A^' lies below (i.e., does not intersect) the cone of influence Qy for every point 
2 G t\ A^. Each such cone of influence corresponds to a causal cone constraint. In ID x Time, a cone of influence is 
defined by two rays with common apex. 

Thus, maximizing the progress of P, and hence the duration of new spacetime elements, requires querying the 
lower hull of cones of influence of points in T \ A^. After the solution is computed on the new front t', including the 
causal slope at every point of t', we obtain a new set of cone constraints. 
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Figure 3.1: A sequence of tent pitching steps in ID x Time. Maximizing the height of each tentpole while staying 
below every cone of influence may require examining remote cones arbitrarily far away. 

For a fixed segment pq incident on p let rsup(p, Q) denote the supremum upper bound on T,+i (p) such that P'Q on 
the front T,+i is causal, i.e., 

T&u^iPiQ) '■= sup{f : P'Q is causal where P' = {p,t)}. 



Let rsup(p) denote the maximum T^xxp{p,Q) for every neighbor Q of P. To maximize the progress at step /+ 1, we 
would like to compute T^upip). The segment P'Q is causal if and only if the slope of P'Q is less than the slope of the 
cone of influence from every point on the front that intersects P'Q. In ID x Time, a cone of influence intersects P'Q if 
and only if the cone intersects the tentpole PP'. 

In general, a cone of influence from arbitrarily far away can intersect the tentpole at p. See Figure [TT] (Nonlocal 
cone constraints do not dominate local cone constraints when the wavespeed everywhere is the same.) Therefore, in 
general, Tsup{p) could be determined by a cone of influence of a point arbitrarily distant from p. In this section, we 
give two algorithms — one to compute Tf^upiPjQ) exactly using ray shooting in an arrangement of rays (lines) and the 
other to approximate rsup(p) numerically using a binary search. 
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Local and nonlocal cone constraints Partition the front into two subsets of points: (i) points in the star 
of P ("local" points), and (ii) points everywhere else on the front ("remote" points). Corresponding to each subset 
we have two disjoint subsets of cones of influence — "ifiocai and "^remote respectively. Each subset of cones limits the 
new time value of p and so the final time value is the smaller of the two values for each of "ifiocai and "demote taken 
separately. 

Consider the subset "^locai- Let CTiocai denote the smallest slope among all cones of influence in "^locai The segment 
P'Q is causal only if its slope is less than Oiocai- Let fiocai be the supremum time value of P' for which the slope of P'Q 



is less than CTiocai- To compute t\oca\ we substitute CTiocai in the condition for causality of P'Q (Equation 2.1 1. 

Next consider the subset "^remote- The front T; is strictly below every cone in 'demote because T; is causal. The 
segment P'Q is causal only if it is also strictly below every cone in "^^remote- Given a cone C S "^remote^ C intersects 
P'Q if and only if C intersects the tentpole PP' . Let fremote denote the smallest time value T for which the tentpole PP' 
where P' = {p, T) intersects exactly one cone in "^remote- The segment P'Q is causal only if T < fremote- 

Therefore, the progress Ti+i{p) — Ti{p) at step /+ 1 is limited because 

Tsu^{p) — min{fiocal,fremote}- 

To maximize the progress at the current step, we choose Ti+\{p) less than rsup(/:>), e.g., equal to T^vi-p{p) minus the 
machine precision. 

We have the following theorem. 

Theorem 3.1. Let X be a causal front and let p be an arbitrary local minimum ofx. Let Wp denote the spatial distance 
from p to its nearest neighbor Then, for every Af such that Q <6X < Wp<J,„i„ the front x' — advance{x , p,At) is causal. 



Proof of Theorem 2.2 Only the segments of the front incident on P = {p, x{p)) advance along with p. Consider an 



arbitrary segment pq incident on p. Since p is a local minimum, we have x{q) > x{p). We have 

x'{p)<x{p)+At 

< x{p)+WpGram 
<x{q) + \pq\a,nm 
<x{q) + \pq\a{P'Q) 

Therefore, the slope of the segment P'Q is less than a{P'Q) and hence P'Q is causal. D 

Computing tREMOTE exactly Computing fremote is equivalent to answering a ray shooting query in the arrange- 
ment of the cones in "^remote- We use a bounding cone hierarchy J^, i.e., a binary tree, obtained from a hierarchical 
decomposition of the space domain to efficiently answer the ray shooting query. The hierarchical decomposition of the 
space domain induces a corresponding hierarchical decomposition of every front. For each element of this hierarchy, 
we store a cone that bounds the cone of influence of every point of the corresponding subset of the front. To answer 
the ray shooting query, we traverse the cone hierarchy from top to bottom starting at the root. At every stage, we store 
a subset "^ of bounding cones such that every cone in "^remote is contained in some cone in the subset '^. The cones 
in 'j^ are stored in a priority queue in non-decreasing order of the time value at which the vertical ray at P intersects 
each cone. Initially, ^ consists solely of the cone at the root of the hierarchy. At every stage, if the cone in '^ that 
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has the eariiest intersection in time does not come from a leaf in the hierarchy then we replace it in the priority queue 
with its children. Continuing in this fashion, we eventually determine the single facet of the front such that the cone 
of influence from some point on this facet is intersected first by the vertical ray at P. The time coordinate of the point 
of intersection is fremote, the answer to the ray shooting query. 

Because the hierarchy is balanced its depth is 0{logm), where m is the number of simplices in the space mesh. The 
tighter the bounding cones, i.e., the better they approximate the actual cones of influence, the better is the efficiency 
of the algorithm. In ID x Time, we observed empirically that on average only a few nodes in the cone hierarchy were 
examined by this algorithm to determine the most constraining cone of influence. 

Approximating Iremote Since we know a range of values [t, (p) + T^i^ , fiocai] that contains fremote, we can approx- 
imate fremote up to any desired numerical accuracy by performing a binary search in this interval. At every iteration, 
we speculatively lift P to the midpoint of the current search interval. Let P" be the speculative top of the tentpole at P. 
We query the cones of influence in "^remote to determine the minimum slope Oiemote among all cones that intersect PP". 
If the maximum slope of the outflow faces incident on P" is less than Oiemote then we can continue searching in the 
top half of the current interval; otherwise, the binary search continues in the bottom half of the current interval. The 
search terminates when the search interval is smaller than our desired additive error A bounding cone hierarchy helps 
in the same manner as before to determine the minimum slope among all cones in ^remote that intersect PP", i.e., all 
cones of influence that contain P". 

3.3 Nonlocal cone constraints in 2D x Time 

In 2D X Time, the algorithm to maximize the progress at each step is complicated by the presence of progress con- 
straints in addition to causality constraints. In this section, we give an algorithm for which the progress guarantee at 
each step is finite and bounded away from zero. The minimum progress guarantee is a function of the local geometry 
of the spatial projection of the front T,- and the global minimum causal slope, analogous to the progress guarantee 
proved in Chapterl2]for the constant-wavespeed case. 

Our algorithm anticipates changing wavespeeds by simulating the next few iterations at each step. The purpose of 
this lookahead is to estimate the actual causal slope encountered in the next few iterations. As a result, we expect that, 
in practice, the actual progress is proportional to the (possibly nonlocal) slope that most constrains causality at the 
current step and in the next few iterations of the algorithm, which may be significantly larger than the global minimum 
slope. Hence, we expect our algorithm to create spacetime elements whose sizes are proportional to the local geometry 
and adapt rapidly to changing causal cone constraints. 

3.3.1 Looking ahead 



From Theorem 2.4 we observe that the new front P'QR is causal if the old front PQR satisfies progress constraint 
o{P'QR). We need to estimate the slope a{P'QR) in the next step in order to enforce a progress constraint in the 
current step. This dependency on the future wavespeed creates the following conundrum. The progress of the front at 
each step / is limited by the progress constraint that must be satisfied by the next front at step / + 1 . However, we do not 
know what the next front is unless we know how much progress is possible at step /. Of course, we can always use the 
global maximum wavespeed as a conservative estimate of the slope a{P'QR); however, this excessively conservative 
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estimate implies a very low progress guarantee for problems where the ratio of maximum to minimum wavespeeds is 
very large. 

To solve this conundrum, we develop a lookahead algorithm that adaptively estimates future wavespeed. For a 
fixed positive integer h, called the horizon, the algorithm at the /th step computes an estimated slope ai,{PQR) for 
every triangle PQR of the current front. The horizon h can be fixed or chosen adaptively at each step. 

When /i = 0, we use the minimum slope Cmin as an estimate of the actual causal slope on the front in the next step, 
so our estimate of future slope is aQ{PQR) = CJmin. When /; > 0, we can use the current estimate to compute the next 
front and the actual slope on this new front to refine our previous estimate. 

Definition 3.2 (/z-progressive). Let h be a nonnegative integer, called the horizon. 

Let PQR be a given triangle. We inductively define APQR as h-progressive as follows. 

Base case h — 0: Triangle PQR is 0-progressive if and only if it is causal and satisfies progress constraint (7,„,„ 
(Definition \2.6^ . 

Case h > 0; Triangle PQR is h-progressive if and only if all the following conditions are satisfied: 

L PQR is causal; 

2. Let P be an arbitrary local minimum vertex of APQR. Let dp denote dist(/7, aff (^r) and let T,„j„ = min{e, 1 — e} 
(Jmindp- Let P'QR = advance (PQR,p, T,„i„) be the triangle obtained by advancing P by T,„i„. Then, PQR must 
satisfy progress constraint u(P'QR) and P'QR must be max{/2 — 1 , 0}-progressive. 

Note that an /i-progressive triangle is also causal. 
Lemma 3.3. For every h > 0, if APQR is h-progressive, then APQR is {h + \)-progressive. 
Proof of Lemma 3.3 If a triangle PQR satisfies progress constraint a (Definition|2.6[), then APQR satisfies progress 



constraint a' for every a' > O. By Lemmas 2.4 and 2.7 if APQR satisfies progress constraint (Tmin, then the triangle 



P'QR after pitching a local minimum P by T^m, where T^m is the quantity in Definition 3.2 is causal and satisfies 



progress constraint CJmin; since Oram is a lower bound on the causal slope anywhere in spacetime, we conclude that 



triangle P'QR is 0-progressive. Therefore, by Definition 3.2 if APQR is 0-progressive, then it is /z-progressive for 



every h>Q. Hence, if APQR is /z-progressive, then it is /i'-progressive for every h' > h. D 

We are now ready to describe our advancing front algorithm in 2D x Time. Recall the causality constraint (Equa- 



tion 



I 2.2 1 and the progress constraint (Definition IZ6]) that applies at each step. 

The causal slope on the new front is computed by the solver; cone constraints can change as a result. In the parallel 
setting, nonlocal cone constraints and updates to the cone hierarchy due to changes in the cone constraints must be 
communicated across processors. 

We claim that being /z-progressive guarantees finite positive progress at each step. If APQR is /z-progressive then 
for every At in the range < Af < min{e, 1 — ejcJrnin width (/?^r) the triangle P'QR obtained by pitching an arbitrary 
local minimum vertex P by Af, is /z-progressive. The amount of progress is a function of Apqr, the distribution of 
wavespeeds, and the parameters e and h. 

Lemma 3.4. Suppose PQR is a triangle of an h-progressive front X for some h > 0. Let P be an arbitrary local 
minimum vertex of APQR. Let dp denote dist(/?, aff^r) and let T,„j„ = min{e, 1 — e\o,nindp. Let P'QR denote the 
corresponding triangle where P' = (p, x{p) + Af ) for an arbitrary Af G [0, Tmm]- 
Then, AP'QR is h-progressive. 
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Input: A triangulated space domain M CV? 

Output: A tetrahedral mesh Q.ofMx [0,°°) 

The initial front Tq is M x {0}, corresponding to time f = everywhere in space. 

Fix h > 0. 

Repeat for / = 0, 1 , 2, . . .: 

1. Advance in time an arbitrary local minimum vertex P = {p,Ti{p)) of the current front T, to P' = 
(p, T,+i(p)) such that T,+i is /i-progressive, and T,+i(/7) is maximized. 

2. Partition the spacetime volume between T and T,+i into a patch of tetrahedra, each sharing the 
tentpole edge PP\ 

3. Call the numerical solver to compute the solution in the spacetime volume between T, and T,+i . 
Figure 3.2: Algorithm in 2DxTuiie with nonlocal cone constraints. 



Proof of LemmaU^ By Definition 3.2 triangle P'QR is {h — 1) -progressive. By Lemma 3.3 triangle P'QR is also 



/z-progressive. D 



By Lemma [34} if the amount of progress made by p in every step is no more than the minimum T^i^ — min{£, 1 — e} 
(Jniindist(/:',aff^r) for every triangle pqr, then every front is /i-progressive. However, the actual progress at each step 
would be no more than the progress guarantee obtained by imposing the progress constraint Omm of Definition 2.6 at 
each step. It is important to minimize the number of spacetime elements by exploiting the fact that the actual wave- 
speed in the next step may be significantly smaller than the global maximum wavespeed; i.e., to take advantage of the 



possibility that o{P'QR) ^ CTmin in Definition 3.2 In the next section, we give an algorithm to greedily maximize the 
progress of local minimum vertex P at each step, possibly at the expense of the progress in subsequent steps, but still 
retain the minimum progress guarantee of T^-am when pitching P. 

3.3.2 Greedily maximizing progress 

We want to maximize the progress at each step in a greedy fashion, i.e., at the /th step given an arbitrary local minimum 
vertex p of the front T, we want to maximize T,+i(p) where T,+i — advance(T,,p, Af) subject to the constraint that Zi+i 
is causal and /^-progressive. 

For a fixed triangle pqr incident on p let Tf,ap{p, QR) denote 

Tsup{p,QR) '■= sup{f : P'QR is causal and /i-progressive, where P' = {p,t)} 

To maximize the progress at step /, we would like to compute T^upip, QR)- 

Similar to the ID x Time case, partition the set of cones of influence from points on the front T, into local and 
remote subsets. Let (Tiocai denote the smallest slope among all local cones of influence. The simplex P'QR is causal 
only if its slope is less than Oiocai- Let fiocai be the supremum time value of P' for which the slope of P'QR is less 



than CTiocai- To compute fiocai we substitute CJiocai in the condition for causality of P'QR (Equation 2.2 1. 

Unlike the ID x Time case, Tsup{p,QR) cannot be computed by ray shooting queries. In 2D x Time, we need an 
oracle to determine which among several cones is intersected for the smallest T by a triangle P'Q'R when the vertex P 
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MaximizeProgress( Front APQR, Vertex P, Integer /; > ): 


1. 


Comment: ah{PQR) is the current estimate of the slope in the next step. 


2. 


a/, <- 


C'lnin i 


3. 


done 


<— false; 


4. 


while not done: 




4.1. 


compute the maximum time value T* such that 

AP'QRwhemP'^{p,T*) 

is causal and satisfies progi-ess constraint ff/,; 




4.2. 


lift p to time value T* giving AP*QR; 




4.3. 


Comment: recursively compute ah-i{P*QR)- 




4.4. 


a' ^ Futures LOPE(P*e/;, h-l); 




4.5. 


done ^- true; 




4.6. 


let P"QR denote the triangle after advancing P 

so that P"QR is causal, it satisfies progress constraint a', 

and the height of PP" is maximized; 




4.7. 


if a{P"QR)> a,,: 

4.7.1. {Comment: Improve the current estimate ff/,.} 

4.7.2. Oh ^ o{P"QR); 

4.7.3. done ^- false; 


5. 


return T*; 



Futures LOPE( Front AABC, Integer h' > ): 


1. 


if/i' = 


= 0: 




1.1. 


return Oraia, 


2. 


CT<— 


X>" 


3. 


for every local minimum, say A, of AABC: 




3.1. 


r' ^ MaximizeProgress(AABC, A, h'); 




3.2. 


let AA'BC be the triangle after advancing A to A' = (a, T'); 




3.3. 


d^min{a,c7(A'BC)}; 


4. 


return ff; 



Figure 3.3: Algorithm to maximize height of tentpole PP' subject to /i-progressive constraints 



of APQR is lifted to f ' — {p, T) while also lifting Qhy a positive amount.. In spatial dimension d >2, the algorithm 
to compute Tsup{p,QR) requires queries involving shooting {d — 1 ) -dimensional faces of the current front. 

Definition |3.2| immediately gives an algorithm to answer the following question: given a triangle PQR in spacetime 
and an integer /z > 0, is APQR /i-progressive? Using the algorithm to answer this decision question, we can approxi- 
mate Tfiup{p,Ql^) up to any given numerical accuracy by performing a binary search in the interval (T,(/7),fiocai] which 
we know contains Tsup{p,QR). The eventual height of the tentpole fP' is at least T^nin, ie., positive and bounded away 
from zero. 



Algorithm MaximizeProgress of Figure 3.3 computes T^upip,QR)- The correctness of the algorithm follows 



from the following observation that FutureSlope estimates the causal slope of the triangle in the next step after 



advancing sufficiently in time. By Lemmas 2.4 and 2.7 FutureSlope advances the local minimum A by at least 
Tmin = min{e,l — e}(7niindist(a,afffec) to A'. Therefore, the estimated slope a computed by FutureSlope is a 
sufficiently accurate estimate of the future slope to ensure that the triangle P'QR obtained by advancing the local 
minimum P to the time value T* computed by MaximizeProgress satisfies the definition of an /i-progressive 



triangle (Definition 3.2 1 
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Maximizing the progress at the current step going from APQR to AP'QR is equivalent to computing the maximum 
slope Gh^P'QR) such that AP'QR is /i-progressive. Algorithm MaximizeProgress relies on a procedure to compute 
the slope anywhere on a triangle A' in the future. We need a procedure to compute the minimum slope among all 
cones of influence emanating from the current front that intersect A'. The algorithm also uses geometric primitives 
such as ray shooting to maximize the height of a tentpole subject to given causality and progress constraints. We will 
postpone discussion of the implementation of these low-level subroutines to a later section. 

We thus have the following theorems. 

Theorem 3.5. For every i > 0, if the front T; is h-progressive, then T; is valid. 



Proof of Theorem \3. 5\ We prove the statement by induction on /. The proof is almost identical to that of Theorem [23] 
except with the additional complication of progress constraints in 2D x Time. The initial front To is progressive by 



definition. By Lemma 2.4 and Lemma 3.4 at each step / the algorithm advances a local minimum vertex p of the the h- 
progressive front T; to the front T,+i such that T,+i is /i-progressive. Let Wp denote the minimum distance dist(/7, aff gr) 
for every edge qr in link(/7). By Lemma 3.4 we know that T,+i > ii{p) + Tmin where T^m = min{£, 1 — £}w pOi-nia- 
Therefore, for every target time T > 0, the entire front achieves or exceeds time T in a finite number of steps. D 



Thus, we obtain a theorem analogous to Theorem 2.10 



Theorem 3.6. Given a triangulation M of a bounded planar space domain where Wmi„ is the minimum width of a 
simplex of M and O is the minimum slope anywhere in M x [0, °°), for every £ such that < £ < 1 our algorithm 



constructs a simplicial mesh of M x [0,7"] consisting of at most 
real T >0, where A is the maximum vertex degree. 



n(T+diam(M)a„, 
min{£, 1 — £}vvv,„>,c 



spacetime elements for every 



Proof of Theorem \3. 6\ By Lemmas |2.4| and |Z7| it follo ws t hat the height of each tentpole constructed by the algorithm 

n(r+diam(M)(j„, 



2.9 



after constructing at most k < 



patches. 



is at least Tmin — min{e, 1 — elvvminCmin- By Theorem : 

the entire front T^ is past the target time T. Since each patch consists of at most A elements, where A is the maximum 

number of simplices in the star of any vertex of M, the theorem follows. D 



3.4 Nonlocal cone constraints in arbitrary dimensions JDxTime 



The algorithm and analysis of Section 3.3 extends in a straightforward manner to higher dimensions d > 3. The 



only additional complications involve the bounding cone hierarchy and queries to this hierarchy. For instance, for 
d — 3, whenever a vertex p of tetrahedron pqrs is advanced in time, we need to query whether the spacetime triangle 
corresponding to Apqr E star(p) intersects any four-dimensional cone of influence. 



We give in this section the definitions in arbitrary dimensions corresponding to those in Section 3.3 and state 
the equivalent theorems without repeating their proofs because they are analogous to the corresponding theorems 
in 2D X Time. 

Consider an arbitrary A:-dimensional face popipi- . .pk incident on pQ. 

Definition 3.7 (/i-progressive). Let h be a nonnegative integer, called the horizon. 

Let PqPiPi- ■ -Pk be a given k-simplex. We inductively define PqP\P2. . .Pk as h-progressive as follows. 

Base case h ~ 0: P0P1P2. . .Pk is 0-progressive if and only if it is causal and satisfies progress constraint Gmin 



(Definition 2.13\. 

Case h > 0; PqPiP2. ■ -Pk is h-progressive if and only if all the following conditions are satisfied: 
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1. PqP\P2. ■ -Pk is causal; 

2. Let Pq be an arbitrary local minimum vertex q/Po^i^2- • Pk- Let dp^ denote dist(/?o, aff /7i/?2- ■ .pk) ond let T,„i„ = 
min{e, 1 — £}(7mindpQ. Let PqPiP2- ■ Pk = advance {PqPiP2. . .Pk,po, T,„j„) be the simplex obtained by advancing 
Pq by Tmin- Then, PqP\P2- ■ Pk must satisfy progress constraint o{Pl^P\P2. . .Pk) and PqPiP2. ■ -Pk must be max{/! — 
1 , 0}-progressive. 

Note that an /i-progressive simplex is also causal. 

Lemma 3.8. Suppose PqPiP2. . .Pk is a simplex of an h-progressive front X for some h > 0. Let Pq be an arbitrary 
local minimum vertex o/PoA^2- • Pk- Let dp^ denote dist(/?0j^ffpi/'2- • -Pk) o-nd let T^i„ ~ min{e, 1 — e}(7,„,„iip(,. Let 
PqPxPi- ■ Pk denote the corresponding simplex where Pq — {po,t{pq) + At) for an arbitrary At G [0, r„„>,]. 
Then, PqP\P2- ■ Pk is h-progressive. 

We have thus shown the finite termination of the algorithm. 

Theorem 3.9. Given a triangulation M of a bounded d-dimensional space domain where Wmin is the minimum width 
of a simplex of M, for every £ such that < £ < 1 our algorithm constructs a simplicial mesh ofM x [0, T] consisting 



of at most 



n{T+diam{M)am, 
min{e, 1 —£\wn^i,^c 



spacetime elements for every real T >0, where A is the maximum vertex degree. 

le height of each tentpole const 
after constructing at most k < 



Proof of Theorem \3. 6\ By Lemmas |2.12| and |2.14[ it follows t hat t he height of each tentpole constructed by the algo 



rithm is at least Tmin = min{e, 1 — e}wminOmin- By Theorem 

patches, the entire front Xk is past the target time T . Since each patch consists of at most A elements, where A is the 

maximum number of simplices in the star of any vertex of M, the theorem follows. D 

3.5 Estimating future slope 

The algorithms of previous sections rely on efficient answers to the following questions: 

1 . Given a triangle PQR and a set of cones of influence, what is the slope of the fastest cone of influence that 
intersects /\PQR1 

2. Given a triangle PQR with P as a lowest vertex and a slope CJ, what is the supremum height of the tentpole PP' 
such that AP'QR has slope less than a? 

Maintaining the entire arrangement of cones of influence is expensive and unnecessary for our purpose; it suffices 
to obtain a cone that bounds (tightly) the actual cone of influence at P. In the absence of focusing, the cone of influence 
of any point P is contained in the cone of influence of every point Q of the front such that P e cone+(2)- 

At each step, we maintain a hierarchical decomposition of the front. We build a bounding cone hierarchy cor- 



responding to this hierarchical partition. See Figure 3.4 We query the cone hierarchy to efficiently maximize the 



progress in time at each tent pitching step. In ID x Time, this query is equivalent to shooting a ray (the infinite ex- 
tension of the tentpole into the future) and determining its earliest intersection with any cone of influence. After each 
step, we update the cones stored in the cone hierarchy to reflect the new wavespeeds computed on the new front. 

Because of no-focusing (Axiom [LT) , we can determine (a lower bound on) the slope at a point P in the future by 
computing all cones of influence from points Q on the current front that contain P. The shallowest such cone deter- 
mines a lower bound on (J(P). It can be computationally very expensive to determine the shallowest cone of influence 
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that contains a given point P. In particular, the shallowest cone of influence containing P may correspond to a point Q 
arbitrarily far from P. To compute this nonlocal cone constraint efficiently, we use a standard hierarchical decomposi- 
tion, called a bounding cone hierarchy, of the space domain. The elements in the hierarchy correspond to subsets of 
the space domain. For each element of the hierarchy, we compute the minimum slope within the corresponding subset 
of the space domain. The smallest element in the hierarchy is a single simplex. In order to determine the strictest 
cone constraint that applies locally, we traverse the hierarchy until we determine the simplex with minimum slope 
whose cone of influence contains P. In practice, we expect that our algorithm has to examine only a small subset of 
the hierarchy; we have observed the resulting speed-up in IDxTime. In the worst case, the algorithm has to examine 
every simplex of the front, but in that case the algorithm will be at most a constant factor slower than one that does 
not use a bounding cone hierarchy. When a patch is solved, the bounding cones are updated with the new slopes by 
traversing a path from a leaf to the root of the hierarchy. This hierarchical approximation technique has been applied 
very successfully to numerous simulation problems, such as the Barnes-Hut divide-and-conquer method (IBH86b for 
A'-body simulations, as well as to collision detection in computer graphics and robot motion planning ( |LMCG96] l and 
for indexing multi-dimensional data in geographic information systems ( IGut84| |. 

Chapter summary 

We have shown how to extend the Tent Pitcher algorithm for arbitrary dimensions to nonlinear problems where the 
wavespeed is not constant. Our expressions for the causality and progress constraints that apply at each step make 
explicit the dependence on the slope of the cone of influence most constraining the progress at that step. This depen- 
dence is not explicit in the formulae of Erickson et al. because they assume without loss of generality that the slope 
is 1 everywhere in spacetime. For the constant wavespeed case, the algorithm in this paper is an alternative to the 
algorithm due to Erickson et al. with potentially weaker progress constraints. We can view the algorithm of Erickson 
et aJ. as looking one step ahead in the sense that the progress constraint at step / guarantees that the front constructed in 
step ! + 1 is causal. Our algorithm can be viewed as looking further — our progress constraint at step / guarantees that 
the front constructed in step / + /z is causal. It needs to be investigated whether the extra complexity of the algorithm 
for h>2 or adaptively choosing h at each step is justified by a more efficient simulation overall. 

The nonlocal nature of the cone constraints pose significant challenges while implementing the algorithm in this 
chapter in parallel. Maintaining a cone hierarchy and performing queries of the sort needed by our algorithm is also a 
significant challenge for spatial dimensions c/ > 3. 
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Figure 3.4: Constructing and traversing a bounding cone hierarchy, (i) All the nonlocal cone constraints limiting the 
height of a tentpole are grouped into (ii) a cone hierarchy (a binary tree) induced by a recursive subdivision of the 
front, (i)-(iv) The cone hierarchy is built bottom-up by merging pairs of bounding cones at each step, (v)-(viii) The 
hierarchy is expanded top-down until the strictest cone constraint is a leaf in the hierarchy. Only a fraction of all cone 
constraints in the hierarchy are examined while traversing the cone hierarchy. 
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Chapter 4 

Adaptive refinement and coarsening 



The duration of spacetime elements constructed by our algorithm is constrained by causality; however, all spacetime 
elements constructed by pitching a simplex of the front have the same spatial diameter as the corresponding simplex 
of the initial front, i.e., the space mesh M. In other words, every front constructed by our algorithm so far is only as 
refined as the initial mesh M. Parts of the spacetime domain, where the solution changes rapidly, need to be meshed 
with smaller elements to achieve the given error tolerance. In this chapter, we make our meshing algorithms adaptive 
to local a posteriori estimates of the numerical error in response to which we refine or coarsen the current front. 
Refinement of the current front means that spacetime elements constructed in future steps will have smaller size, and, 
after a finite number of refinement steps, reduced error. By adapting the size of the spacetime elements to the error 
estimate, we are able to compute a more efficient mesh for a given error bound. Without adaptivity, we would require 
the initial space mesh to be as refined as the finest resolution required at any time, which would necessitate many more 
elements for meshing the same spacetime volume. 

In this chapter, we concentrate on the case of constant slope (i.e., constant wavespeed) where CT denotes the 
slope everywhere in spacetime or more generally a — (Jmin is the reciprocal of the globally maximum wavespeed. In 
ChapterlSlwe will discuss adaptive refinement and coarsening in the presence of changing wavespeeds. 

Also, in this chapter, we rectify an oversight in the proof of correctness in our paper (lACE^04l l by including a case 
in the proof of Theorem |4. 3 [ that was missing from the corresponding theorem in the paper. 

4.1 Problem statement 

Mesh adaptivity in our advancing front framework is the problem of adapting the spatial size of front facets. We 
require an algorithm that responds to error tolerance requirements by reducing the spatial size of front facets where 
necessary so that future spacetime elements are smaller and therefore have reduced error When the error estimate is 
sufficiently greater than the allowed error, we want the algorithm to increase the spatial size of front facets wherever 
permitted so that future spacetime elements are bigger Reducing the error is a requirement for the solution strategy to 
advance; hence, it is necessary for the algorithm to respond immediately to refinement demands. However, coarsening 
is a request — an excessively refined mesh still gives a sufficiently accurate solution but it is not efficient. The algorithm 
should coarsen aggressively to be efficient but it can decide to postpone a coarsening request until a later step. 

Our main challenge is to incorporate refinement and coarsening of the front at each step into our advancing front 
framework. Each iteration of the advancing front algorithm, i.e., each application of the advance function, is under- 
stood to be a tent pitching step followed by zero or more refinement and coarsening operations. The actual operations 
for refining and coarsening the current front are described in detail later in this chapter 

Mesh refinement and coarsening affects spacetime elements in the future when some vertex of the currently refined 
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front is pitched. Since our algorithm coarsens only a pair of triangles that have been previously obtained by refinement, 
it is advantageous to start with the coarsest acceptable space mesh, e.g., the coarsest mesh that is still a good enough 
piecewise linear approximation to the space domain, and let the meshing algorithm adapt the mesh locally as necessary 
in response to the demands for numerical accuracy. In Chapter [7] we discuss generalization of coarsening beyond a 
simple undoing of a previous refinement. 

In dimension d >2, we saw in Chapter l2] that nontrivial progress constraints are necessary to guarantee positive 
progress at each step. These progress constraints are functions of the shape of the triangles on the front. But the 
shape of front facets, or more accurately the shape of the spatial projection of the front facets, is changing as a result 
of refinement and coarsening! Therefore, the challenge is to modify the progress constraint to anticipate changes in 
shape due to an arbitrary amount of refinement. 

Our solution For one-dimensional space domains, we have proved in Chapter l2] that every (causal) front is 
valid. Therefore, modifying Tent Pitcher to make it adaptive in ID x Time is very easy. The input space mesh is 
a one-dimensional simplicial complex. To refine the front, we bisect a front segment into two parts, say two equal 
halves; coarsening reverses this operation by merging two segments that are coUinear in spacetime. Since refinement 
and coarsening does not alter the gradient of the time function restricted to the elements involved in these operations, 
refinement and coarsening preserve causality. 

Note that when pitching a local minimum vertex P, an interior vertex, the smaller of the temporal aspect ratios 
of the resulting two triangles is maximized if the two edges incident on P have approximately equal lengths of their 
spatial projections. We use this fact to guide our choice of where to subdivide a front segment. 

In higher dimensions, we define progressive fronts and prove that if a front is progressive then it is valid. For 
li = 2, we give an algorithm that, given any progressive front T,, constructs a next front T,+i — advance (T,-,p, A?) such 
that T,+i is progressive and T,+i {p) is maximized. Whenever p isa local minimum of T;, the progress T,+i (p) — Ti{p) 
is guaranteed to be at least Ty^m, which is a function of the input and bounded away from zero. It is necessary to predict 
the shape of the spatial projections of triangles on the new front T,+i and ensure that T, satisfies progress constraints 
that anticipate refinement and coarsening of the new front. We choose a refinement method, called newest vertex 
bisection, that allows us to predict all possible shapes of triangles on any front in the future after an arbitrary number 
of refinement and coarsening steps. We incorporate constraints associated with every such shape in the progress 
constraints that must be satisfied by the front at each step, and call such a front a progressive front. The details and 
formal definitions are in Section !?. 2.41 

Our algorithm adapts to a posteriori estimates of numerical error as follows. A patch is solved as soon as it is 
created. If the estimated numerical error, i.e., estimated energy dissipation, for any element in the patch is greater 
than some threshold ^i then the element is marked for refinement and the patch is rejected. If no element is marked 
for refinement, then the patch is accepted. If the largest estimated numerical error for an element is less than some 
threshold £,2, where ^2 < ^b then the element is marked as coarsenable. If a patch is rejected, then the front is not 
advanced and for every element marked for refinement, the corresponding facet of the current front is bisected. Since 
repeated bisection decreases the spatial diameter, the size of future spacetime elements decreases; we assume that a 
finite number of refinement steps eventually reduces the numerical error. 

For now, we consider only the case d = 2, i.e., space domains M that are 2D simplicial complexes consisting of 
triangles, edges, and vertices embedded in some ambient Euclidean space. Adaptive refinement and coarsening in 



dimensions d > 2 remains an open problem and we postpone a discussion of higher dimensions until Section 4.3 
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Input: A one-dimensional space mesh M C E' 

Output: A triangular mesh H of M x [0,oo) 

The initial front To is M x {0}, corresponding to time f = everywhere in space. 

Repeat for / = 0, 1 , 2, . . .: 

1. Advance in time an arbitrary local minimum vertex P = {p,Ti{p)) of the current front T,- to P' = 

(p,Ti+i(p)) such that T,+ i is causal and Ti+i{p) is maximized. 

If all segments of T; incident on P are marked as coarsenable by the solution during the 
previous iteration, then choose P' so that the coarsenable segments become collinear in 
time, unless the height T,+ i(/?) — T,(/7) of the resulting tentpole is too small. 

2. Partition the spacetime volume between T and T,+i into a patch of triangles, each sharing the tentpole 
edgePP'. 

3. Compute the solution in the spacetime volume between T, and T,+i as well as the a posteriori error 
estimate. 

4. If the spacetime error indicator tolerates the error in the patch, then advance the front to T,+i and 
merge every adjacent pair of coarsenable collinear segments. 

5. Otherwise, some segments of the front are marked for refinement; bisect each such segment to get 
the new front T,+i. 

Figure 4.1: Adaptive meshing algorithm in ID x Time 

If a patch is rejected, then for every tetrahedron PP'Q'R' marked for refinement the triangle pqr in the spatial 
projection is bisected. A pair of triangles can be coarsened by merging them if the result is a single triangle. If the two 
triangles were previously obtained by bisecting a larger triangle, then they can always be merged whenever they are 
coplanar in spacetime. 



Figure |4T| describes our adaptive meshing algorithm in ID x Time; figure [435] describes the adaptive algorithm 
in 2D X Time. Both algorithms proceed by pitching local minima, and they refine and coarsen the front in response 
to error estimates. In the parallel setting, we repeatedly choose an independent set of local minima of the current 
front, equal to the number of processors, to be advanced in time simultaneously. The resulting patches can be solved 
independently. If a patch is accepted, the local neighborhood of the front is advanced without any conflicts with other 
patches. Thus, the solution everywhere in spacetime is computed patch-by-patch in an order consistent with the partial 
order of dependence of patches. 

4.2 Meshing in 2D x Time 

In this section, we give new adaptive progress constraints in 2D x Time that guarantee that each front, even after an 
arbitrary number of refinement and coarsening steps, is able to make positive progress by tent pitching. 
First, we describe our mesh refinement procedure. 
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4.2.1 Hierarchical front refinement 

Mesh refinement is needed to achieve the desired numerical accuracy. Intuitively, the larger is the variation in the 
solution within a subdomain, the smaller is the desired element size. Where the solution is changing slowly and 
smoothly, it is sufficient to use larger elements — since the number of elements is smaller, the solution procedure is 
more efficient for the given numerical accuracy. Therefore, the mesh should be only as refined as necessary. 

For tracking evolving phenomena efficiently, it is also necessary to coarsen previously refined elements when the 
phenomenon recedes and the solution is no longer changing rapidly. Thus, coarsening is desirable for computational 
efficiency. However, over-refinement does not hurt the quality of the solution. 

Since our DG method uses a discontinuous formulation, we do not require a smooth grading of element size. In 
fact, experimental results suggest that if the mesh successfully tracks a sharp shock, then it is perfectly acceptable to 
use very coarse elements aligned with the shock trajectory. 

In the terminology of the different kinds of adaptivity that have been studied in the literature — p-adaptivity (adapt 
the degree of the polynomial basis functions), /z-adaptivity (adapt the size and number of elements in the mesh), and 
r-adaptivity (adapt the locations of nodes in the mesh, e.g., smoothing) — we perform /z-adaptivity to adapt the size of 
spacetime elements. 

The choice of bisection method There is a vast body of Hterature on mesh refinement for both structured 
and unstructured meshes. See the survey by Jones and Plassman ( JP97) for an introduction. Specifically, we are 
interested in hierarchical mesh refinement by recursive bisection. Regular subdivisions of a simplex were studied by 
Bank et aJ. (IBSW83I I in 2D, and Edelsbrunner and Grayson dEGOOI l in higher dimensions. Mesh refinement in 2D 
requires bisecting specified triangles to decrease their diameter as well as propagating to neighboring triangles to 
maintain a triangulation. The longest edge refinement introduced by Rosenberg and Stenger (IRS75I I has been later 
popularized by Rivara ( Riv97l lRiv84l |Riv96l l. Newest vertex bisection was originally developed by Sewell (ISew72l) 
and later adapted by Mitchell jMit88t IMit89t IMit91l) . Arnold et al. (.AMPOQj , Bey (Bey95l ), Biinsch jBan91l ). Liu 
and Joe JLJ94I ILJ95I) . Maubach JMau95b . and Bey (BeyOOl extended newest vertex bisection to higher dimensions. 



Most importantly, Arnold et al. jAMPOOI l proved that the number of shapes generated by their recursive bisection of a 
simplex is finite. 

Our adaptive algorithm uses the newest vertex bisection refinement method, originally developed by Sewell (ISew72b . 
later adapted by Mitchell (IMit88l IMit89l IMit91l) in the context of multigrid methods, and still later studied and gen- 
eralized to three dimensions by Bansch ( |Ban9U . This method is similar to, but not identical to, longest-edge refine- 
ment lRS75llRiv84llRiv96l) . 

We call the newest vertex of a triangle its apex and the opposite edge its base. Initially, one vertex of each triangle 
in the mesh is chosen arbitrarily as its apex. Newest vertex bisection replaces a triangle with two smaller triangles, 
each with half the area of the original triangle, obtained by bisecting along the line segment through the apex and the 
midpoint of the base. The new vertex introduced at the midpoint is the newest vertex of both smaller triangles. See 



Figure 4.2 



The descendants of any marked triangle under newest vertex bisection fall into only eight homothetic classes 



(Figure 4.2 1. There are only four directions along which the triangle or any of its descendants can be further bisected. 
These four directions are parallel to the three edges of the triangle and the bisecting segment. Each of the four ways 
to choose three of the four directions gives two similar triangle shapes that are mirror images of each other, for a total 
of eight triangle shapes, equivalent up to translation and scaling. Any two triangles in the same equivalence class have 
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Figure 4.2: Newest vertex refinement ( lACE+041 . Newest vertices of each triangle are marked. 

corresponding newest vertices. Refinement by three levels is guaranteed to decrease the diameter by at least a half. 

A front triangle PQR and its grandparent AABC create homothetic tetrahedra when their corresponding vertices 
are pitched. This is because causality and progress constraints are gradient constraints that scale with the spatial 
geometry of front triangles. Therefore, our refinement algorithm using newest vertex bisection adapts the resolution 
of the spacetime mesh but is limited in its ability to adapt the shape and temporal aspect ratio of spacetime tetrahedra. 

Propagation When we refine one triangle in a mesh, we may be forced to refine other nearby triangles in order to 
maintain a conforming triangulation. Maintaining a triangulated front is necessary for our algorithm — nonconforming 
or dangling vertices create complications that we do not know how to solve yet. 

We call an edge of the triangulation a terminal edge if it is the base of every triangle incident on it. 



Suppose vertex a is the apex of a l\abc that is bisected (Figure 4.3 i. If be is not a boundary edge then some 
neighboring triangle IScbe shares the edge be. To maintain a triangulation, ISebe must be bisected also. If be is not 
the base of Aebe, then the child of Aebe sharing edge be must be bisected as well, and the bisection of Aebe will 



propagate recursively; see Figure 4.3 



It is easy to prove that this propagation terminates, regardless of which vertex is chosen as the apex of each 
triangle. Suppose a refinement of Aabe with apex a propagates to a neighboring triangle sharing the edge be. A single 
newest vertex bisection of Aabe makes the edges ab and ae the bases respectively of the two children of Aabe. If the 
propagation revisits Aabe, it must return by bisecting either the edge ab or the edge ae. But since these edges are now 
terminal edges, the propagation terminates at that step. 

The propagation path touches every triangle in the worst case, but in practice, the propagation path usually has 
small constant length; see Suarez et al. (ISPC03I I for an analysis of a similar refinement algorithm. Because we bisect 
triangles first and then propagate, instead of the other way round, our algorithm terminates, regardless of which vertices 
are marked in each triangle: each triangle in the mesh is bisected at most twice. (Mitchell's original head-recursive 
algorithm (IMit88t |Mit89l |Mit91| l can enter an infinite loop if the initial selection of marks does not satisfy a certain 
matching condition.) 

The propagation follows a directed path in the dual graph of the triangulation, leaving each triangle along the dual 
edge corresponding to the base of the triangle. The propagation terminates when the base of the last triangle bisected 
is a terminal edge. The propagation either encounters a boundary edge or the propagation revisits one of the children 



of a triangle that was bisected earlier along the propagation path; see Figure 4.4 

Newest vertex bisection never subdivides any angle of a triangle more than once. If the apex of each triangle in 
the initial mesh is the vertex with largest angle, then every vertex v is the apex of at most 5 triangles incident on v. 
(In the degenerate case where six equilateral triangles meet at a vertex, we can break ties using a straightforward 
symbolic perturbation scheme.) Thus, the degree of a vertex v in the initial space mesh increases by at most 5 as a 
result of refinement. If v is not a vertex of the initial space mesh then the degree of v is 4 when v is inserted and is at 
most 8 at any subsequent step. We conclude that if a triangulation has maximum degree A, then any refinement of that 
triangulation has maximum degree at most max{A + 5,8}. 
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(a) (b) 

Figure 4.3: Refining ISabc propagates to neighboring triangles: (a) before and (b) after refinement. 




(13) (14) 

Figure 4.4: Refinement propagation path terminates when it revisits a triangle. 
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Figure 4.5: Bisecting an edge and the two triangles incident on it 

BlSECTTRlANGLE(Triangle Apqr): 

1 . Let qr be the base of Apqr, 

2. Let Aqrs be the neighbor of Apqr sharing the edge qr, 

3. Comment: If gr is a boundary edge, there is no such neighbor. 

4. Bisect the edge qr by inserting a new vertex m at its midpoint; 

5. Insert the edges pm and ms; 

6. Mark all four angles at m; 

7. If qr is not the base of Aqrs: 

8. Mark Aqrs; 



Figure 4.6: Bisecting a triangle and its neighbor. 

Basic mesh operations The entire refinement propagation can be expressed as a sequence of edge bisection 
and edge flip operations. Bisecting a terminal edge achieves newest-vertex bisection of all triangles incident on it. 
However, if an edge e is not a terminal edge, then bisecting e leaves at least one of the triangles incident on e in a 
dirty state because an edge of this triangle other than its base has been bisected. To rectify this situation, we clean the 
dirty triangle by another edge bisection followed by an edge flip. In this fashion, every clean triangle is subdivided 
according to the rules of newest-vertex bisection. Dirty triangles are transient and are cleaned before the surrounding 
neighborhood of the front is advanced. 



Edge bisection This operation bisects an edge shared by one or two triangles; see Figure 4.5 The new vertex is the 



apex of all four new triangles. Bisecting a triangle at level I produces two new children at level l+l; every triangle in 
the original space mesh is at level zero. The set of all triangles that ever appear on any front form a forest of rooted 
binary trees with each triangle in the initial space mesh as the root of some tree in the forest. The height of a node in 
the refinement tree is the number of edges on a longest path from the node to any of its descendants. 
Edge flip An edge flip replaces one diagonal of a convex quadrilateral with the other diagonal. Note how in Fig- 
ure [47] the apex of the adjacent triangle (shown dashed) is changed. On the left, we see Apqr first split into Apqu 
and Auqr, then Auqr is split into Auqv and Auvr. On the right, after the flip operation, we have Apqr first split 
into Apqv and Apvr, and then Apvr split into Apvu and Auvr. Vertex p is the apex of Apqr. The refinement subtree 
rooted at Apqr changes as a result of the edge flip. 
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Figure 4.7: Edge flip 



Lazy propagation Edge bisections caused by propagation unnecessarily create smaller triangles in parts of the 
front even where the error estimate is sufficiently below the acceptable threshold. If the propagation were delayed until 
a later step, the spacetime mesh would contain fewer tetrahedra. In a parallel setting, interprocessor communication is 
inefficient; so, it is useful to allow the algorithm to proceed on unrelated processors and not require all processors to 
wait until the entire propagation sequence has been executed. 

Therefore, we propagate refinements lazily, by temporarily bisecting an edge of a neighboring triangle other than its 
base and marking this triangle as dirty. A dirty triangle can be cleaned up later by performing one edge bisection, which 
in turn can be propagated lazily, plus one edge flip. Lazy propagation is particularly useful for parallel implementation, 
where the mesh may be distributed over multiple processors, because it minimizes both communication costs and 
deadlocks. If a processor needs to refine a triangle in its subset of the mesh, it need not wait for the entire propagation 
sequence (which may be controlled by other processors) to be completed before proceeding with its next task. Even 
in the serial case, lazy refinement may reduce the number of changes to the mesh, since later coarsening may stop the 
propagation of refinement early. 
Earnest vs. lazy propagation Algorithm Refine is called to refine a triangle currently on the front. By defini- 



tion, this triangle is a leaf in the refinement forest. The contrast between earnest propagation (Figure 4.8 1 and lazy 



propagation (Figure 4.9 1 is shown in Figure 4.12 Earnest propagation causes the adjacent triangle to be split which 



may propagate further. Lazy propagation (bottom path) splits the adjacent triangle but does not propagate further 



immediately — the gray arrows in Figure 4. 12 are not traversed until later Instead, the propagation is delayed until the 
adjacent triangle needs to be refined or pitched. In the interim, the triangulation may consist of transient dirty triangles. 
After cleaning up, the final result is the same. 



Refinement witli lazy propagation 

As long as all dirty triangles are cleaned up, refinement with lazy propagation is identical to refinement with earnest 



propagation; see Figure 4.12 
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4.6 1 



REFlNE(Triangle Apqr): 

1. Refine ANDPROPAGATE(Apg'r); 

REFlNEANDPROPAGATE(Triangle Apqr): 

1. Let e be the base of A/^gr; 

2. Let T be the neighboring triangle th at sh ares the edge e; 

3. BlSECTTRlANGLE(A/:>^r); (Figure ' 

4. If T is marked: 
4.L CleanUp(t); 

CLEANUP(Triangle Apqr): 

1. Let /7 be the apex of Apgr; 

2. Let Apqs and Asqr be the current children of Apqr; 

3. {Comment: The edge opposite q is currently bisected.} 

4. REFINEANDPROPAGATE(Ai^r); 

5. Flip the edge qs; 



Figure 4.8: Refinement with earnest propagation. 



4.2.2 Coarsening 



Coarsening is the opposite of refinement and hence called de-refinement; we coarsen locally by undoing a single edge 
bisection. Unlike refinement, coarsening does not require propagation further into the mesh to maintain a conforming 
triangulation, although one coarsening step may make other coarsenings possible. In particular, if we refine a triangle 
and then immediately coarsen, we can (but need not) coarsen along the entire refinement path. 

DeRefine is called to coarsen a triangle currently on the front that is a strict descendant of some triangle in the 
original mesh. By definition, the triangle being coarsened is a leaf in the refinement forest and its parent exists. See 



Figure 4.13 



When the propagation path of a refinement loops back on itself, there will be no degree-4 vertex in the portion of 



the mesh corresponding to this loop (Figure 4.14i, and algorithm DeRefine will not apply. The solution is to perform 
an edge flip to create a degree-4 vertex and then call DeRefine. The edge flip will create a dirty triangle, and this 
dirty triangle can be coarsened immediately by DeRefine along with its neighbor. 

4.2.3 Adaptive meshing in 2D x Time 

We are now ready to describe an iteration of our advancing front algorithm. Advance a single vertex p by a positive 
amount, where p is any local minimum of the current front T,-, to get the new front T,+i such that every triangle pqr 
on Ti+i is progressive. In the parallel setting, advance any independent set of local minima forward in time, each 
subject to the above constraint. 

To incorporate adaptivity into our meshing algorithm, we make a small change to the main loop. At each iteration, 
just as before, we choose a local minimum vertex P of the front, move it forward in time to create a tent, and pass the 
tent and its inflow data to the spacetime DG solver. We assume that the solver also computes an a posteriori estimate of 
its own numerical error If the error within any element of the tent is above some threshold, the solver rejects the patch, 
at which point the meshing algorithm throws away the tent and refines the facets of the front whose elements had high 
error. (Alternately, we could refine the facets until their diameter is smaller than a target length scale computed by the 
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LAZYREFlNE(Triangle Apqr): 

1. If parent(A/?^r) exists and is marked dirty: 

1.1. {Comment: Apqr is transient and so we cannot refine Apqr.} 

1.2. CLEANUPl(parent(A75g'r)); 

2. Else: 

2. 1. LAZYREFINEANDPROPAGATE(A/?^r); 

LAZYREFlNEANDPROPAGATE(Triangle Apqr): 

1 . Let e be the base of Apqr, 

2. Let T be the neighboring triangle th at sh ares the edge e; 

3. BlSECTTRlANGLE(A/:>^r); (Figure ' 



4.6 1 



4. If parent(T) exists and is marked dirty: 

4.1. {Comment: T is transient.} 

4.2. CLEANUP2(parent(T)); 
Cleanup 1 (Triangle Apqr): (Figure 4.10 1 

1 . {Comment: Apqr is marked and has heigiit 1 .} 

2. Let p be the apex of Apqr; 

3. Let Apqs and Asqr be the current children of Apqr, 

4. {Comment: The edge opposite q is currently bisected.} 

5. LAZYREFINEANDPROPAGATE(A5q'r); 

6. Flip the edge qs; 
CLEANUP2(Triangle Apqr): (Figure 4.11 1 

1. {Comment: Apqr is marked and has height 2.} 

2. Let p be the apex of Apqr, 

3. Let Apqs and Asqr be the current children of Apqr, 

4. {Comment: The edge opposite q is currently bisected.} 

5. If Apqs is a leaf: 

5.1. {Comment: Asqr is subdivided.} 

5.2. Flip the edge qs; 

6. Else if Asqr is a leaf: 

6.1. {Comment: Apqs is subdivided.} 

6.2. Let Apst and Atsq be the children of Apqs; 

6.3. REFINEANDPROPAGATE(Ai^r); 

6.4. Flip the edge qs; 

6.5. Flip the edge st; 

Figure 4.9: Refinement with lazy propagation. 

solver) Note that this refinement may propagate far outside the neighborhood of P. We accept the numerical solution 
and update the front only if the error within every element of the patch is acceptable. 

On the other hand, the error estimate within an element may fall below some second threshold, indicating that 
the mesh is finer than necessary to compute the desired result. In this case, the DG solver marks the outflow face of 
that element as coarsenable. We can coarsen four facets of the front into two only if they are the result of an earlier 
refinement, they are all marked as coarsenable, and each pair of triangles to be merged is coplanar To make coarsening 
possible, our algorithm tries to make coarsenable sibUngs coplanar, by lowering the top of the tent. However, to avoid 
very thin elements, we accept the lower tent only if its height is above some threshold. If the lower tent is accepted 
and its outflow faces are still marked coarsenable, then we coarsen the front. This means, of course, that merging a 
pair of coarsenable triangles may be delayed due to the coplanarity constraint imposing a too-small tentpole height. 
We definitely observe this in practice, leading to a more inefficient mesh than absolutely necessary. In fact, coarsening 



64 




Figure 4.10: Before (solid) and after (dotted) one edge bisection and an edge flip performed by CleanUpI 

P P 





(a) 



(b) 



Figure 4.11: Before (solid) and after (dotted) CleanUp2: (a) only an edge flip is required; (b) an edge bisection 
followed by two edge flips are required 

may be delayed several steps. In Chapter [6] we give an algorithm to ensure that every coarsenable pair of triangles 
will eventually be made coplanar and merged, while still guaranteeing a bounded minimum tentpole height. 



See Figure 4.15 for an outline of the adaptive Tent Pitcher algorithm. See Figure 4.16 for an illustration of the 



effect of refine and coarsen operations interspersed with tent pitching steps. 

Cleaning up after a lazy propagation If refinement is being propagated lazily, we must ensure that the 
triangles incident on p are clean. 

CleanupBeforePitching( Vertex p): 

1 . For every triangle /\pqr incident on p\ 

2. If parent (Ap^r) exists and is marked dirty: 

3. Cleanup l(parent(Ap^r)); 



COPLANARITY CONSTRAINTS While pitching p, we see whether coarsenable triangles adjacent to p can be made 
coplanar with their siblings after pitching. For each triangle /\pqs incident on p, consider the four triangles neighbor- 
ing /\pqs referred to in DeRefine. If all the four triangles are currently marked as coarsenable, we solve for x'{p) to 
compute the new value of T'(p) that would make p coplanar with q, r, and u. This is the value of t'{p) that satisfies 
what we call the coplanarity constraint. Satisfying the coplanarity constraint may require a shorter tentpole at p\ if 
the shorter height is less than some constant (such as i/io) times the height without the constraint, then we choose to 
ignore the coplanarity constraint to avoid creating degenerate elements. Otherwise, the final height of the tentpole at p 
is hmited by the coplanarity constraint and coarsenable triangles can be merged immediately after pitching p. 
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Figure 4.12: Earnest propagation (top path) vs. lazy propagation (bottom path) 

4.2.4 Causality and Progress Constraints 

To complete the description of our adaptive version of Tent Pitcher, all that remains is to define progressive triangles 
and fronts, i.e., to define the progress constraints that limit the new time value for each vertex to be pitched. Each 
triangle in the space mesh imposes constraints on the time values at each of its vertices. When we pitch a tent over a 
local minimum vertex P, the new time value T,+i(/5) is simply the largest value that satisfies the constraints for every 
triangle pqr incident on p in the space mesh. 



One constraint is of course the causality constraint of Equation 2.2 repeated below: 



\PP\ 



< 



Vct^-IIVt/ 



cir\ 



The more subtle and important change in our algorithm is the introduction of new progress constraints. In the 
earlier non-adaptive algorithm (IEGSU02| |. the progress constraint was a function of the shape of the underlying space 
elements. In our new algorithm, the shape of the underlying element is subject to change; each triangle in the current 
front may be the result of any number of refinements since the last time any ancestor or descendant of that triangle 
was pitched. Consequently, our progress constraints must take into account the shapes of all possible descendants of a 
triangle simultaneously. This requirement motivates our choice of newest vertex refinement; because the descendants 
of any triangle fall into only a finite number of homothety classes, we have only a finite number of simultaneous 
progress constraints. 

We can visualize these constraints by referring back to our circle diagram; see Figure |4. 18| A minimal set of 
progress constraints can be obtained by simply taking the union of the non-adaptive progress constraints of the eight 
shapes in the hierarchy. Recall that each of the eight shapes, corresponding to a descendant of Apqr, defines a 
forbidden zone where the gradient vector of Apqr must not lie. We call these the "primary" forbidden zones. The 
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DEREFlNE(Triangle Apqs): 

1. Let Apsr be the sibling of Apqs, i.e., edge qr is currently bisected; 

2. Verify that Apqs and Apsr are leaves and are coplanar; 

3. Verify that s has degree 4; 

4. Let u be the fourth neighbor of s other than {p, q, r}; 

5. Verify that Arsu and Ausq are leaves and are coplanar; 

6. {Comment: Equivalently, r-s-q must be coUinear.} 

7. {Comment: All four triangles incident on s need NOT be coplanar.} 

8. Replace Apqs and Apsr with their parent Apqr; 

9. Replace Arsu and Ausq with their parent Arqu; 

10. {Comment: The two new triangles Apqr and Arqu are not dirty.} 
n . If parent( Ap^r) exists and is coarsenable: 

12. DEREFWEiApqr); 

13. If parent(Ar^M) exists and is coarsenable: 

14. DEREFlNE(Ar^M); 



Figure 4.13: Algorithm to de-refine 




Figure 4.14: De-refining a loop in the refinement propagation path requires edge flips to create degree-4 vertices. To 
merge the first pair of triangles, an edge flip is necessary, (1)^(2), to create a degree-4 vertex that can now be deleted. 
Subsequent coarsening steps (2)^(3) and (3)^(4) each require an edge flip. 
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Input: A triangulation M C M 

Output: A tetrahedral mesh of M x [0,oo) 

Initial front Tq is the space mesh at time f = 0, i.e., To{p) = for every p G V{M) 

Repeat for / = 0, 1 , 2, . . .: 

• Lift a local minimum vertex p of the current front T, to obtain the new front Tj+i such that every 
triangle pqr on T,+i is progressive. 

• Solve the resulting patch. 

• If the patch is rejected by solver, one or more elements in the patch are marked for refinement. 

• For every element marked for refinement, bisect its inflow facet; otherwise, advance the front. 

• If any pair of coarsenable siblings are coplanar, then coarsen this pair. 



Figure 4.15: Adaptive Tent Pitcher algorithm with refinement and coarsening 



(a) 




Figure 4.16: A sequence of 5 tents pitched by the adaptive algorithm: (a)-^(b) pitch twice, (b)— >(c) refine and pitch, 
(c)^(d) pitch, (d)^(e) pitch, (e)^(f) coarsen and pitch 

progress constraint that the corresponding triangle PQR on the front must satisfy is the conjunction of the progress 
constraint for each of these eight shapes. In addition, we must also forbid "secondary zones" from which the only 
alternative, due to the nature of Tent Pitcher, is to enter a primary forbidden zone. Staying outside the primary and 
secondary forbidden zones is sufficient to guarantee that the progress at each step is bounded away from zero. 

Fix two real numbers £ and (p such that 0<e<^<(l + e)/2<l. For any triangle Aabc with apex a, we define 
the diminished width of Aabc as follows. 



Definition 4.1 (Diminished width (lACE"'"04t '). Let abc be an arbitrary triangle, where a is the apex. (See Figure 4.17.) 
The diminished width of Aabc is defined as 

(1 — e)dist(fl,aff/7c), 

v{abc) :— min ^ (1 - ^)dist(/7, aff ac), 

(1 — ^)dist(c, affafe) 



The first distance is measured from the apex to the opposite edge and is scaled differently from the other two 
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Figure 4.17: Aabc where a is the apex. 



Z^ ^ -^ 





Figure 4.18: The adaptive progress constraint is obtained from the progress constraint for each of the eight different 
shapes of descendant triangles. 

altitudes. This definition extends recursively to any descendants of Aabc obtained by newest vertex subdivision; in 
the interest of readability, we will always list the vertices of any triangle with the apex first. We express our new 
progress constraints algebraically by limiting the difference in time values along each edge in the subdivided triangle, 
as follows. 

Definition 4.2 (Adaptive progress constraint a). Fix e and (p such that 0<e<^<(l + e)/2<l. Let AABC be 
an arbitrary triangle of a front T with apex A. Let d, e, and f be midpoints of sides be, ac, and ab respectively 



(Figure 4.17). We say that the triangle ABC satisfies adaptive progress constraint a if and only if 



\T{a)-T{b)\ <2w{dca)(J 
|T(fl) -t(J)| < w{abc)a 
|T(fl)-T(c)| <2w{dab)a 
\x{b)-x{c)\ <4w{ead)a 



The constraints in the definition of the progress constraint apply recursively to all descendants of the triangle Aabc, 
but these recursive constraints are equivalent to one of the four constraints above. 

Now suppose we are pitching a triangle Apqr of the front T, where T{p) < T{q) < T(r). Let w be the midpoint 



of qr; let u be the midpoint of pr; and let v be the midpoint of pq; see Figures |4.19| and |4.20| Depending on which of 
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the three vertices is marked as the apex, the new time value t'(/7) is bounded as a result of these progress constraints 
in the three different ways enumerated below. Notice that when p is not the apex, lifting p also lifts either u or y, so 
progress constraints along edges qu or rv also indirectly limit l' {p). 



If /? is the apex: 



T{q) +2w{rpw)a^ 
T'(/5)<min<^ 'z{s) + w{pqr)a, 
'z{r) + 2w{wpq)o 



(4.1) 



If r is the apex: 



T{q)+Aw{uvr)a, 
t'{p) <mini^ 'z{r) + 2w{vqr)a, 

2x{r) - %{q)+2w{rpq)a ^ 



(4.2) 



If ^ is the apex: 



t'{p) < min 




2w{uqr)G, 
4w(wuq)(y, 
T(r) +2w{rpq)a 



(4.3) 



The adaptive progress constraint (Definition |4.2|l is stronger than the non-adaptive progress constraint (Defini- 



tion 



2.6 1; in particular, the new progress constraint limits the gradient of all three edges of APQR, not just the highest 



edges, and anticipates future refinement and coarsening of triangles on the front by newest vertex bisection by also 
limiting the gradient along the bisector edge. Of course, our algorithm works even when no refinement or coarsening 
operations are actually performed. 

The following theorems were proved by Abedi et al.. In this section, we rectify an oversight in the proof of 
correctness in the paper by Abedi et al. ( lACE+04l l by including a case in the proof of Theorem|4.3|that was missing. 



Theorem 4.3. If a front t,- is progressive, then for every local minimum vertex p, for every Af e [0, EWpOmin], the front 
T,+i = advance(Ti,p,At) is causal. 

Proof. Since only the triangles of the front incident on P advance along with p, we can restrict our attention to an 
arbitrary triangle pqr incident on p. Let T and t' denote T,|„„,- and T,+i| . respectively. Let p be the orthogonal 
projection of p onto line qr. 

We consider two cases separately. 



Case 1: t(p) > T(q) > T(p) See Figure 4.19 In this case, we have Zpqr < '^ji. We will show that 

|lVT|^J|<(l-e)c7 
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Figure 4.19: Triangle pqr with /.pqr < ^/i with apex at p, q, and r respectively. 



i.e., that 



T(r)-T(^)<(l-e)|^r|c7. 



The proof then follows from Lemma 2.4 because the antecedent of the lemma is satisfied. 



We consider three cases depending on which vertex of Apqr is the apex. 
If p is the apex, then 



T(r) — t(^) < 4w{upw)G 

< 4(1 -e)dist(M,aff/5w)a 
= 2{l - e)dist{r,aff pw)a 

<2(l-e)|vvr|c7 
= {l-e)\qr\a 



If r is the apex, then 



x{r) - -l{q) < 2w{vrp)a 

<2(l-e)dist(v,aff/?r)(7 
= (1 — e)dist(g',affpr)(7 

<(l-e)H(j 



If q is the apex, then 



T(r) — x{q) < 2w{upq)a 

< 2(1 — e)dist(M, aff/?g')(7 
= (1 — e)dist(r, affp^)a 

<(l-e)H(j 



(by the progress constraint) 
(by definition of diminished width) 

(2dist(M,affpw) = dist(r, aff/7w)) 

(dist(r, affpw) < \wr\) 

(2\wr\ = \qr\) 



(by the progress constraint) 

(by definition of diminished width) 

(2dist(v, aff /?r) =dist{q,aSpr)) 

(dist(g', aff pr) < \qr\) 



(by the progress constraint) 

(by definition of diminished width) 

(2dist(M,aff/7(7) = dist(r, aff/jg')) 

(dist(r, affpg') < \qr\) 
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(i) (ii) (iii) 

Figure 4.20: Triangle pqr with j^pqr > ^/i with apex at p, q, and r respectively. 



Case 2: t(p) < T(q) See Figure 4.20 In this case, we have Zpqr > '^/i. We will show that 

||VT|^J<(l-e)cTsinZMr 

i.e., that 

T{r)-T{q) < (l-£)|^r|(JsinZ/7^r. 

The proof then follows from LemmalZ4]because the antecedent of the lemma is satisfied. 



See Figure 4.20 a). Let j3 — \pq\/\pp\- Since \pq\ ^ 0, we have 

T'ip)-T{p) _T'{p)-T{q) , T{q) - Tip) \pq\ 



\PP\ 



\PP\ 
Z'jp) - Tjq) 

\PP\ 



\pq\ \pp\ 



Using Equation 4.4 the causality constraint (Equation 2.2 1 can be rewritten as 

T'ip) - X{q) 



\PP\ 



<j0^-\\^A,r\? 



-iSIIVTl II 



We consider three cases depending on which vertex of /\pqr is the apex. 
If p is the apex, then 



x{r) — x{q) < Aw{upw)o 

<4(1 — e)dist(M,aff/?w)a 
= 2(1 -e)dist(r,affpw)o■ 



(by the progress constraint) 

(by definition of diminished width) 

(2dist(M,aff/7w) = dist(r, aff/7w)) 



(4.4) 



(4.5) 
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Hence, 



, dist(r,affnw) 

\wr\ 

= (1 —e)asinZpwr 
\PP\ 



(l-e)(7 



<(l-e)a 



|/?w| 
\PP\ 



\Pl\ 
= (1 — e)G sin Zpqr 



because \pw\ > \pq\ 



If r is the apex, then 



T(r) — x{q) < 2w{vrp)a 

<2(l-e)dist(v,aff/?r)(7 
= (1 -e)dist(^,aff/?r)o' 



(by the progress constraint) 
(by definition of diminished width) 

(2dist(v, aff /5r) = dist(^,affj3r)) 



Hence, 



z{r) — T{q) dist{q,affpr) 



m 



^ \^ 


^'^ \^r\ 


= (1- 


-e)asmZqrp 


= (1- 


\pr\ 


<(1- 


\pq\ 


= (1- 


-e)a sin Zpqr 



because \pr\ > \pq\ 



If q is the apex, then 



x{r) — x{q) < 2w{upq)a 

< 2(1 — e)dist(M,aff/?^)c7 
= (1 — e)dist(r, aff/?^)^ 



(by the progress constraint) 

(by definition of diminished width) 

(2dist(M,aff/5^) = dist{r,affpq)) 



Hence, 



T(r) - T(g) ^ dist(r,aff/7g) ^ 

\qr\ - \qr\ 

= (1 — e)G sin{n — Zpqr) 
= (1 — £)a sin Zpqr 
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n 

Theorem 4.4. If a front T; is progressive and if < £ < (p < {I + £) /2 < 1, then for every local minimum vertex p 
there exists T,„j„ > 0, where T,„i„ is a function of the triangle Apqr and the parameters £ and (p, such that the front 
T,+i = advance {Ti,p, At) is progressive for every At G [0,7;,„„]. 

Proof. Since only the triangles of the front incident on P advance along with p, we can restrict our attention to an 
arbitrary triangle pqr incident on p. Let T and t' denote T,| . and T,+i| . respectively. Let p be the orthogonal 
projection of p onto line qr. 



Consider the progress constraints (Equations 4.1-4.3 1. Each progress constraint limits the progress T'{p) — T{p) 



from above. The only progress constraint for which this limit is not obviously positive is Equation |4. 3 1 which applies 
when q is the apex of Apqr, i.e., 

t:'{p) < 2x{q) - x{r)+2w{qrp)a 

Subtracting i^p) from both sides, because x{q) > i{p), the above constraint is satisfied if 

T(r) — t(^) < 2w{qrp)o 

Since APQR is progressive, we have 

x{r) — x{q) < 2w{upq)o 

Therefore, it suffices to show that w{qrp) — w{upq) is positive and bounded away from zero. Expanding the definition 
of diminished width, we rewrite X = w{qrp) — w{upq) as X = min{A,B, C} where 

{(1 — e)dist(M,affpq'), 
(1 - (p)dist{p,aff qu), 
(1 -(p)dist{q,aff up) 

> (1 — e)dist(^,aff/7r) — (1 - (p)dist(^,affM;7) 
= {(p - £)dist{q,aS pr); 

{(1 — £)dist{u,aSpq), 
(1 -^)dist(;?,aff^M), 
(1 — (p)dist{q,aff up) 

> (1 - (p)dist{r,aff pq) — {I - £)dist{u,affpq) 
= 2(1 - (p)dist{u,aff pq) - (1 — e)dist(M,aff/?q') 
= (1 +e — 2^)dist(M,aff/?q'); 

{(1 — e)dist(M,aff/?^), 
(1 -(p)dist{p,affqu), 
(1 -9)dist(^, affM/?) 

> (1 — 9)(dist(/7,aff^r) — min{dist(/:i,affq'M),dist(^,affM/?)}). 
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Since £ < ^ < (l + e)/2, we have A > andB > 0. See Figure |4.1 9 iii) or Figure 4.20 lii); we also have 



dist(p,aff,.)^?5'^i;^ 

2area{Aupq) 
\uv\ 
2arsa{Aupq) 



> 



max{|M/?|,|M^|} 



The inequahty above follows because the bisector segment uv must be shorter than at least one of the two sides up 
and uq. Now, 2area(AMp^)/|Mp| = dist(^,affpr) = dist(g', aff wp) and 2area(AM/?^)/|M^| = dist(p,aff^M). Hence, it 
follows that C> 0. n 

4.2.5 Changing the apex 

Sometimes it may be advantageous to change the apex (newest vertex) of a triangle. For instance, if we wish to bisect 
a given edge without causing the refinement to propagate to neighboring triangles, we can make the edge a terminal 
edge by choosing the apexes of incident triangles appropriately. In Chapter IT] when we discuss pitching inclined 
tentpoles to account for boundary motion or to perform smoothing, the apex of a triangle may no longer correspond to 
the largest angle of the triangle. It is useful to subdivide the largest angle, or equivalently to bisect the longest edge, 
to guarantee sufficient progress in the next step. We therefore require some mechanism to change the apex or newest 
vertex of a given triangle. 

The adaptive progress constraint can be strengthened to allow for changing the apex of any triangle. Specifically, 
for each triangle pqr, we impose all three types of progress constraints on Apqr when p is the apex, when q is the 
apex, and when r is the apex simultaneously. Thus, we say that a causal triangle ABC is progressive if and it satisfies 



the adaptive progress constraint of Definition 4.6 for every choice of apex of the triangle. 

Hence, when pitching a vertex p of Apqr we impose all the following constraints simultaneously on the new time 
value t'(p): 

T{q)+2w{rpw)(7, t{s) + w{pqr)a, T{r) + 2w{wpq)a, 
t:'{p) <mm^ T{q)+4-w{uvr)a, T{r) + 2w{vqr)a, 2T{r) - T{q)+2w{rpq)a, ^ (4.6) 

T{q)+2w{uqr)a, T{r) +4-w{wuq)a, z{q) — T{r) +2w{rpq)a 



Since the new progress constraint is stronger. Lemma 2.4 still holds 



Lemma 4.5. If a front T is progressive and if < £ < (p < {I + e) /2 < I, then for any local minimum vertex p 
there exists A > 0, where A is a function of the triangle Apqr and the parameters £ and (p, such that the front 
t' — advance(x,p,M) is progressive for every Af G [0, A]. 



Proof. The proof is almost identical to that of Lemma 4.5 which does not rely on the fact that the three cases. 



depending on which vertex of triangle pqr is the apex, are mutually exclusive. D 
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4.2.6 Another look at the adaptive progress constraint 

The progress constraint of Abedi et al. is equivalent to the following constraint for every descendant triangle by newest 
vertex bisection. Abusing notation slightly, we refer to an arbitrary descendant as Aabc with apex a as in Figure |4TT7l 
Then, the bisector ad of Aabc must have gradient bounded by 

|T(a) - T:{d)\ < w{abc)G. 

Expanding the definition of diminished width and dividing both sides by \ad\, we obtain an equivalent progress con- 
straint: 

\T{a) — T{d)\ , ( dist(fl,aff/7c) dist(fe,affflc) dist{c,affab)\ 



len(aJ) [ len(flii) ' len(flii) ' len(flii) J 

Observe the following identities: dist(^,affflc) = 2dist(ii, aff ac); dist(c,affafe) = 2dist{d,afiab); and sinZadb = 
sin /Lade. Hence, we can rewrite the last inequality as 

IJV f|adll < niin{(l -e)sinZa(ic, 2(1 -9) sin Zfifac, 2(1 - (p)s\nAdab}o 

Therefore, the adaptive progress constraint <7 due to Abedi et al. that limits the gradient of Aabc along each of the 
four edges in the subdivided triangle can be rewritten equivalently as follows: 

Triangle ABC satisfies progress constraint a if and only if each of the following conditions is satisfied: 

1. II V t|„^|| = ||V t|j^|| < min {(1 — e) sinZdea, 2(1 — 9) sin Zeda, 2(1 — 9) sinZeiic}c7 

2. II V t|^^/|| < min {(1 — e) sin Zadc, 2(1 — 9) sinZdac, 2{l~(p) sin Zdab} a 

3. ||Vt|„J| = ||VT|^,y.|| < min {{I - e) sinZd fb, 2{l - (p) sinZf db, 2{l - (p) sinZf da} a 

4. ||V z\J\ = ||V t|^J < min {(1 -e) sinZfga,2{l - cp) sinZgfd,2{\ - <p) sinZgfa}a 



Observing similarity and supplementarity of angles in Figure |4.17| we can rewrite the above constraints to refer 
only to angles in Aabc and its two children. 

Definition 4.6 (Restatement of adaptive progress constraint). Triangle ABC with apex a satisfies the adaptive progress 
constraint if and only if each of the following conditions is satisfied: 

II V AubW ^ ™" {(1 ~£) sinZbac, 2{l — (p) sin Zbad, 2{l~(p) sin Zabc} a 
II V '^ladW - ™i" {(1 ^£) sinZadc, 2{l — (p) sinZdac, 2{l — (p) sin Zdab} a 
W^ '^lacW - ™i" {(1 ^^) sinZbac, 2(1 — (p) sinZacb, 2{l — (p) sin Ziiacjff 
II V t|,^^.|| < min {(1 -e) sin ZMa, 2(1 — 9) sinZacb, 2(1 — 9) sin Zabc} a 



Lemma 4.7. Triangle APQR satisfies tlie progress constraint if and only if its every ancestor and every descendant 
also satisfies the progress constraint. 

Proof. Since (p < (l + e)/2, we have 2(1 — (p) > (1 — e). Therefore, because of the correspondence between the angles 



of Apqr and those of an ancestor or a descendant, Aabc satisfies the progress constraint (Definition 4.6 1 if and only 



if every ancestor and every descendant of Apqr by newest vertex bisection also satisfies the progress constraint. D 

76 



4.3 Mesh adaptivity in higher dimensions 

Generalizations of newest vertex bisection to higher dimensions exist and are well-known. For instance, the refinement 
method of Maubach (IMau95l IMau96b generates a bounded number of similarity classes of simplices, independent 
of the amount of refinement due to repeated bisection. The main challenge we confront is to extend our progress 
constraints in 3D x Time to account for adaptive refinement and coarsening. We anticipate more complications than 
the interplay of primary and secondary forbidden zones which we encountered in 2D x Time. 

Our progress constraint still allows us to perform limited smoothing and retriangulation of the front at each step. 

Under admittedly very strong constraints on the gradient of the front, it is possible to adapt the front to improve 
the mesh resolution and temporal aspect ratio of future spacetime elements. 

Lemma 4.8. If every angle 9 in the spatial projection of a front x satisfies sin0 > <^„n„, andif\\V x\\ < (1 —£)<j),„i„(7, 



nun*~^m[n> 



then any front t' obtained by retriangulating T is causal and satisfies progress constraint Omin of Definition 2.13 



Proof. Since CTmin is a lower bound on the reciprocal of the fastest wavespeed anywhere in spacetime, T is clearly 
causal. 

Since ||Vt|| < (1 — e)0min<7min; for any spatial direction vector n we have (Vt) -n < ||Vt|| < (1 — e)0minO'mm- 
Therefore, for any face F in the spatial projection of t', we obtain || V t|^|| < (1 — e)^!!!^^!!!^- D 

Thus, when the gradient of the front T is bounded, essentially any arbitrary front modification is permitted as long 
as the weak simplicial complex property is maintained; in other words, a retriangulation T -^ t' is permitted only if t' 
conforms to T or if the operation can be achieved by adding new non-degenerate tetrahedra each of which has a causal 
outflow facet. 

Chapter summary 

Adaptive meshing is crucial for efficient simulations. Unless the size of mesh elements adapts to changing require- 
ments, an extremely fine mesh would be required throughout the domain. The number of elements required to be 
solved in such a fine mesh would render the problem intractable. In this chapter, we presented an adaptive meshing 
algorithm that incorporates refinement and coarsening into our advancing front framework. Our algorithm adjusts, as 
soon as required, the size of future spacetime elements to a posteriori estimates of the energy dissipation after solving 
each patch. 

For the purposes of the meshing algorithm, we interpret the dissipation-based error indicator as a black box, 
specifically, a binary valued function pL that given a spacetime element E either accepts or rejects E. To allow the 
algorithm to coarsen the mesh, the lower and upper thresholds on either side of the target dissipation should be 
sufficiently well-separated by a hysteresis zone. 

Our algorithm strives for an acceptable balance between accuracy and efficiency. Aggressive coarsening, for the 
purpose of minimizing the number of spacetime elements, must be balanced with the minimum acceptable element 
quality needed for accurate numerical analysis. 

We refer the reader to the paper by Abedi et aJ. ( lACE+041 for a sample application to a scientific problem of 
practical complexity: the phenomenon of crack-tip wave scattering within an elastic solid subjected to shock loading. 
As expected, our experimental results show a significant improvement in the quality of the solution as a result of 
adaptivity, especially near discontinuities in the solution or its derivative. Also, we were able to achieve a better 
solution without using a fine mesh everywhere, which would have resulted in a massive increase in computation time. 
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We would like to extend all the current results to higher dimensions, specifically to 3D x Time. For the actual mesh 
refinement, we can directly use the generalization of newest vertex bisection to three dimensions by Bansch (|Ban91| ) 
and others. The main theoretical hurdle in higher dimensions is deriving minimal progress constraints that guarantee 
that the front can always advance. Another concern is to be able to describe the constraints in such a fashion that they 
are easy to implement. It seems likely that analytic constraints should suffice in higher dimensions, so the algorithm 
to maximize progress at each step can be solved exactly. 
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Chapter 5 

Adaptive Meshing with Nonlocal Cone 
Constraints 



In Chapter|4] we gave a meshing algorithm for 2D x Time that refines and coarsens the front to adapt the size of future 
spacetime tetrahedra in response to numerical error estimates. The algorithm of Chapter p] assumed a fixed worst- 
case estimate of the causal slope at each step. In Chapter l3] we relaxed this conservative estimate. The result was 
an algorithm that adapted the duration of spacetime tetrahedra to changing wavespeeds. In this chapter, we give an 
algorithm to simultaneously adapt the mesh resolution to error estimates and anticipate changing wavespeeds due to 
nonlinear response. 

We now consider the problem of supporting adaptive refinement and coarsening in 2D x Time via newest vertex 
bisection in the presence of changing wavespeeds. As noted in Chapter|4] devising necessary and sufficient progress 
constraints to support refinement and coarsening via higher-dimensional analogues of newest vertex bisection remains 
an open problem. The progress constraint at the /th step depends on the geometry of the front, which may change in 
step /+ 1 due to refinement and coarsening. The amount of progress going from T, to T,+i is constrained by the fastest 
wavespeed on the new front T,+i. As before, we would like to characterize the front T, as valid if and only if it is 
guaranteed to make finite progress to T,+ i such that the new front T,+i is also valid. The problem is to anticipate a fast 
wavespeed in the future simultaneously with changes in the front due to refinement and coarsening. Our contribution 
is to combine the progress constraints of Chapters [2] and [4] with the lookahead algorithm of Chapter |3] We define a 
new notion of progressive fronts and prove that every progressive front is valid. 

The front is refined in response to numerical error estimates. There is another reason for the meshing algorithm 
to refine the front. A spacetime element created by pitching a vertex of a front T may encounter vastly different 
wavespeeds. Note that the corresponding triangle of T must satisfy a progress constraint that depends on the wavespeed 
encountered after pitching. Hence, if the maximum wavespeed on any outflow facet of the element is much larger than 
the maximum wavespeed on any inflow facet, then the temporal aspect ratio of the tetrahedron is small. The height of 
the element is proportional to the ratio of maximum to minimum wavespeed anywhere in the element. If the temporal 
aspect ratio is too small, the element may be unsuitable for computing the numerical solution. A spacetime element 
with very small temporal aspect ratio may have a singular stiffness matrix, which may make it impossible to solve, 
or may require too many iterations for the solution to converge. If the solution does converge, the element will likely 
accumulate a large error and will therefore be rejected; the cost of solving the element will be incurred without any 
advancement of the front. There is empirical evidence that the accuracy of the solution within an element is positively 
correlated with the quality of the element, defined as the ratio of its inradius to its circumradius; a small decrease in 
the height of the spacetime element may imply a much larger decrease in its quality. 

The computational cost of solving an element is significantly higher than the amortized per-element cost of mesh 
generation. DG solvers for nonlinear PDEs may take several iterations for the approximate solution to converge. If the 
accuracy of the solution and its efficiency can be improved by a better quality mesh with fewer elements, the additional 
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cost of mesh generation is expected to be easily compensated by the reduced total computation time. 

Therefore, for correctness as well as efficiency, it is important for the meshing algorithm to avoid creating space- 
time elements with small temporal aspect ratio whenever possible. 

In this chapter, we give a meshing algorithm that refines the front in response to a large increase in causal wave- 
speed, independent of and in addition to any error-driven refinement demanded by the solver. While simulating the 
next few tent pitching steps to estimate future wavespeed, if the algorithm predicts a bad quality element due to in- 
creasing wavespeed, then it can preemptively refine the front in the current step. Refining a triangle improves or 
maintains up to a constant factor the temporal aspect ratio of the tetrahedron created by pitching the triangle. Also, 
more progress can be made in the current step if the algorithm prepares for future refinement 

In practice, especially for examples where the Mach factor — the ratio of global maximum to minimum wavespeeds— 
is not too large, it may be acceptable to allow the solver alone to guide front refinement. We cannot yet claim an 
improvement in temporal aspect ratio for arbitrary wavespeed distributions, as a result of the algorithm in this chapter. 
Regardless, we believe that this algorithm is of theoretical interest. The algorithm solves a two-parameter greedy 
optimization problem and illustrates the tradeoff between the two parameters — one a horizon parameter that measures 
future wavespeeds and the other an adaptive lookahead parameter that estimates the need to refine the front in the 
future. Thus, the meshing algorithm presented in this chapter attempts, in practice, to reduce the number of tetra- 
hedra and improve their worst-case quality. Experiments are needed to determine whether the complexity of a more 
complicated algorithm is outweighed by a better quality and/or more efficient solution. 

Recall from Chapter [3] that the algorithm maintains a bounding cone hierarchy to speed up the computation of the 
most restrictive cone constraint at each step. Whenever the front is refined or coarsened, we update the hierarchical 
decomposition by recomputing the hierarchy for a subset of the front. For the sake of efficiency, we rebalance the 
hierarchy after a significant refinement or coarsening of any subset of the front. Instead of updating the slope and 
number of cones in the hierarchy immediately after each step, we can update in a lazy fashion because the old cone 
constraints are still valid conservative estimates of new constraints. Lazy updates are likely to be efficient in a parallel 
setting by reducing interprocessor communication. 

5.1 Adaptive meshing in ID x Time with nonlocal cone constraints 

We note that it is trivial to perform adaptive refinement and coarsening of the front in ID x Time in the presence of 
nonlocal (possibly anisotropic) cone constraints. We maximize the height of each tentpole subject only to the constraint 
that each front is causal, as described in ChapterlS] 



The proof of the following theorem is almost identical to that of Theorem 2.2 



Theorem 5.1. Let x be a causal front and let p be an arbitrary local minimum of X. Then, for every Af such that 
< Af < Wp(7,ni„, the front t' = advance (T,p, At) is causal. 



Proof of Theorem 5.1 Only the segments of the front incident on P = (p, 'z{p)) advance along with p. Consider an 



arbitrary segment pq incident on p. Since p is a local minimum, we have t(^) > T(p). We have 

x'{p)<T{p)+M 

< T(p) + Wp(Jmin 

<Tiq) + il-d)\pq\a{P'Q) 
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Hence, 

II l,K/ii i^^i V ^^ 

Therefore, P'Q is causal. D 

Note that when pitching a local minimum vertex P, an interior vertex, the smaller of the temporal aspect ratios of 
the resulting two spacetime triangles is maximized if the two edges incident on P have comparable spatial lengths. 
Therefore, our adaptive strategy is guided by the desire to keep the spatial lengths of adjacent front segments roughly 
equal. 

Because of Theorem |5.1[ we obtain a lower bound on the worst case height of any tentpole. Thus, we have proved 



Theorem 2.3 in the situation where the wavespeed is not constant throughout spacetime. 



5.2 Adapting the progress constraint to the level of refinement in 2D x Time 

In this section, we consider the adaptive meshing problem in 2D x Time in the presence of changing wavespeeds. 
We restrict ourselves to mesh adaptivity by newest vertex bisection of front triangles and the reverse operation. See 
Chapterfflfor a discussion of newest vertex bisection. 

The key property of newest vertex bisection is that each descendant of a triangle PQR has its edges parallel to three 
of the four directions defined by APQR and its apex — the three edges PQ, QR, and PR plus the bisector PS through 
the midpoint s of the base qr. 



Recall from the statement of Lemmas 2.4 and 2.7 that the front at each step must satisfy a progress constraint that 
depends on the wavespeed in the next step. When the wavespeed is not constant, Gmin may not be a good estimate 
of the future wavespeed. We have improved the situation somewhat by giving an algorithm to compute an estimate 
^h{PQR) of the wavespeed encountered in the next step by APQR of the current front, by looking ahead h steps for 



some horizon parameter h. Thus, we can use the estimated slope ah{PQR) as the radius of the disc in Figure 2.2 when 
imposing adaptive progress constraints on APQR. However, this choice of ah{PQR) may be too conservative because 
it uses the same wavespeed O), {PQR) in the progress constraint on APQR as on each of its descendants after bisection. 
A key property we use in deriving the new algorithm in this chapter is that the boundary of a cone is a ruled surface. 
Recall from Chapter |3] that we partition the cone constraints while pitching A into local and nonlocal constraints. 
Therefore, if any triangle AABC intersects a given nonlocal cone, then this intersection is a line segment in the plane 
of AABC. Therefore, if the bisector edge AD intersects a cone of influence, then so do at least two of the edges of 



AABC (Figure 4.17 1. Therefore, if a fast wavespeed in the future intersects AD but does not intersect both edges AB 
and AC of AABC, then it can do so only if AABC has been bisected at least once. 

Thus, our contribution is to make the progress constraint imposed on a given triangle due to the shape of its 
descendants adaptive to the level of refinement of the descendant triangles. 

The progress constraint that we impose on a given triangle PQR of a front T depends on two parameters, h and I. 



The parameter h is the horizon, as in Definition 3.2 For each descendant AABC of APQR by / newest vertex 
bisections, we enforce a progress constraint on APQR that depends on the shape of AABC, the level of refinement I, 
and the horizon parameter h. The greater the level of refinement of a triangle A, the greater is likely to be the estimated 
slope Gf,{A) for a fixed h, i.e., if A2 is a descendant of Ai by newest vertex bisection, then 0/,(A2) > a/j(Ai). 
Therefore, we anticipate up to / levels of refinement where / > is another parameter. For a triangle A 1 at refinement 
level greater than I, we impose a progress constraint that depends on its ancestor A2 at level /; in other words, we 

81 



draw a figure similar to Figure 2.2 except with a disc of radius d'/,(A2). If / = 0, then A2 = APQR; the advantage of 
choosing / > is that in practice we expect a/,(A2) ^ di,{PQR) because A2 is smaller than APQR. In other words, 
for the same horizon parameter h, the smaller triangle A2 will encounter a subset of the cone constraints in the next h 
steps of pitching vertices of A2 as the larger triangle APQR. For the same number of lookahead steps h, a smaller 
triangle is likely to advance through more tent pitching steps than its larger ancestor because each tent pitching step 
makes progress in time proportional to the size of the triangle. 

To support adaptive refinement and coarsening by newest vertex bisection, we have to enforce the adaptive progress 
constraint of Definition |4.2| We therefore define adoptively h-progressive fronts (Definition |3.2| l, to refer to this 
stronger adaptive progress constraint, as follows. 

Definition 5.2 (Adaptively /i-progressive). Let h be a nonnegative integer, called the horizon. 

Let PQR be a given triangle. We inductively define APQR as adaptively /j-progressive as follows. 

Base case h = 0: Triangle PQR is adaptively 0-progressive if and only if it is causal and satisfies the adaptive 
progress constraint Omin (Definition\4.2\. 



Case h > 0: Triangle PQR is adaptively h-progressive if and only if all the following conditions are satisfied: 

L PQR is causal; 

2. Let P be an arbitrary local minimum vertex of APQR. Let dp denote dist(;?, aff gr) and let T,„i„ = min{£, 1 — e} 
Oniindp. Let P'QR — advance{PQR, p, T,„j„) be the triangle obtained by advancing P by T„ji„. Then, PQR must 
satisfy the adaptive progress constraint o{P' QR) and P'QR must be maxj/i — 1 , 0}-progressive. 

We say that a front T is adaptively h-progressive if every facet of T is adaptively /i-progressive (Definition|5.2[i. 



The following two lemmas are analogous to Lemmas 3.3 and 3.4 



Lemma 5.3. For every h > 0, if APQR is adaptively h-progressive, then APQR is adaptively {h+ l)-progressive. 
Proof of Lemma\5J\ If a triangle PQR satisfies adaptive progress constraint a (Definition |4.2[i, then APQR satisfies 



progress constraint a' for every a' > (7. By Theorems 4.3 and 4.4 if APQR satisfies the adaptive progress con 



straint ffmin, then there exists Tmm > such that the triangle P'QR after pitching a local minimum P by Tmin is causal 
and satisfies the adaptive progress constraint dmin- Since Omm is a lower bound on the causal slope anywhere in space- 



time, we conclude that triangle P'QR is adaptively 0-progressive. Therefore, by Definition 5.2 if APQR is adaptively 
0-progressive, then it is adaptively /j-progressive for every h>Q. Hence, if APQR is adaptively /z-progressive, then it 
is adaptively /I'-progressive for every h' > h. D 

Lemma 5.4. Suppose PQR is a triangle of an adaptively h-progressive front X for some h > 0. Let P be an arbitrary 
local minimum vertex of APQR. Let dp denote dist(/7, aff^r). 

Then, there exists r„„„ > such that AP'QR where P' = [p, x{p) + Af)/or an arbitrary Af G [0, r,„,„] is adaptively 
h-progressive. 



Proof of Lemma^^ By Definition 5.2 triangle P'QR is adaptively [h — 1) -progressive. By Lemma 5.3 triangle P'QR 



is also adaptively /j-progressive. D 

Definition 5.5 ((/z,/)-progressive). We inductively define a triangle APQR as (h,l)-progressive if it is h-progressive 



(Definition 3.2) and each of the two children obtained by newest vertex bisection of APQR is (/i,max{/ — 1,0})- 
progressive. 



Base case I = 0: APQR is {h,0)-progressive if it is adaptively h-progressive (Definition 5.2 \. 
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We say that a front is {h,l)-progressive if every triangle on the front is (/!,/)-progressive. We say that a front is 
progressive if it is (/z, Z)-progressive for some fixed h and /. 

A progressive front is guaranteed to make sufficient positive progress after at most a finite number of refinements 
by newest vertex bisection of its facets. 

Lemma 5.6. If a front T is (h,Q)-progressive, then for every local minimum p of x there exists a T,„i„ > 0, such that 
the front t' — advance(t,p,At) is (h,0)-progressive for every At € [0, r„„„]. 



Proof of Lemma p^ By definition, every triangle PQR of the front T is adaptively /i-progressive. By Lemma 5.4 the 
triangle P^QR obtained by advancing local minimum vertex P to P' by Tmin is also adaptively /z-progressive. Hence, 
there exists a T^ia > such that the front t' = advance(T,/:',Af) is (/i,0)-progressive for arbitrary Af G [0,rn,in]. □ 



By Definition 5.5 if a triangle PQR is (/i,/)-progressive, then any triangle obtained by a single newest vertex 
bisection of PQR is (/!,/') -progressive for /' = max{/— 1,0}. 

We relax the definition of a valid front (Definition |2T| to permit refinement of the front one or more times at each 
step without advancing the front. 

Definition 5.7 (Vahd front). We say that a front T is valid if there exists a positive real At bounded away from zero such 
that for every T G 5R- there exists a finite sequence of fronts T, Ti, T2, . . ., Tk where T^. > T, each new front T,+i in the 
sequence obtained from the previous front T,- by one of two operations: either (i) advancing some local neighborhood 
of Ti by At, or (ii) newest vertex bisection of some triangle ofTj. 

Therefore, we conclude the following theorem. 

Tlieorem 5.8. If a triangle PQR is {h,l)-progressive, then APQR is valid. 

Proof of Theorem \5. 8\ The proof is by a very simple induction on the parameter /. If I = 0, then the theorem follows 



from Lemma 5.6 Otherwise, if / > 0, then by the induction hypothesis, the child triangle ABC of APQR, which is 



(/z,/ — 1) -progressive, is valid. Hence, by Definition 5.7 the theorem follows for the case I > also. D 



5.2.1 A unified adaptive algorithm in 2D x Time 

Our unified algorithm greedily maximizes the progress such that each front is (/z,/) -progressive for some choice of h 



and /. The algorithm can be as complicated as desired. Definition 5.5 stresses the fact that our algorithm can optimize 
the choice of h and I, likely doing better than the theoretical guarantee obtained by setting h = l = 0; however, a simple 
choice of h = I = I may perform sufficiently well in practice. 

5.2.2 A simpler algorithm in 2D x Time 

Arguably, the algorithm presented in the previous sections is very complicated to implement for general h and I. Is this 
complication really necessary? In the case of discontinuities in wavespeeds due to shocks, we think the complication 
is unavoidable. However, in the case that the wavespeed function is smooth, we present a simpler algorithm that is 
guaranteed to terminate. For each triangle A on the front, estimate the maximum and minimum wavespeed at any point 
of A. If the ratio of the maximum to minimum wavespeeds is too large, say greater than 4, then refine the triangle; 
otherwise, if this ratio is too small, say less than 2, then mark the triangle as coarsenable. Since the wavespeed is 
continuous, this refinement will stop after a bounded number of steps. If the wavespeed function is discontinuous, we 
would need to impose a lower bound on the element size to limit refinement. 
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Input: A triangulated space domain M C E^ 

Output: A tetrahedral mesh of M x [0,oo) 

Initial front To is the space mesh at time f = 
Fix h, I. 

Repeat for / = 0, 1 , 2, . . .: 

1. Advance in time an arbitrary local minimum vertex P = {p,'^i{p)) of the current front T, to P' = 

{p, T,+i (/?)) such that T,+i is (/z, Z)-progressive, and T,+i (/?) is maximized. 

2. Partition the spacetime volume between T and T,+ i into a patch of tetrahedra, each sharing the 
tentpole edge PP\ 

3. Solve the resulting patch. 

4. If the patch is accepted, then some outflow triangles on the new front T,+i are marked coarsenable. 
Coarsen any pair of coarsenable sibling triangles if they are coplanar on the new front T,+i, as long 
as the resulting coarser front is also (/z,/)-progressive. 

5. If the patch is rejected, then one or more triangles on the front T, are marked for refinement. Per- 
form newest vertex bisection of every triangle marked for refinement, propagating the bisection to 
neighboring triangles to maintain a triangulated front. 

Figure 5.1: Unified adaptive algorithm in 2D x Time. 

Chapter summary 

Refinement of a triangle on the front is easy to incorporate into the advancing front meshing algorithm with nonlocal 
cone constraints because any faster cone of influence that intersects a smaller triangle also intersects its parent and is 
accounted for in the progress constraint that must be satisfied by the parent. 

Coarsening presents a challenge when the wavespeed is not constant everywhere. When two smaller triangles are 
coarsened, the maximum wavespeed on the larger triangle is the maximum of the wavespeeds of the original smaller 
triangles. Therefore, both triangles being coarsened must satisfy the same progress constraint before they can be 
coarsened, i.e., we can merge two triangles into one if and only if the new triangle is /z-progressive. When coarsening 
is possible only under such constraints, we need to carefully prioritize each coarsening step so that the front is only as 
refined as necessary and not much more. Prioritizing the various coarsening requests remains an open problem. 



84 



Chapter 6 

Target time conformity 



For comparing the accuracy of the numerical resuhs with those obtained by other techniques, such as for convergence 
studies, we sometimes need a spacetime mesh of the bounded volume from time f = to f = T for some target 
time r > 0. Also, some modifications of the front, such as coarsening by derefinement, are permitted only when the 
gradient of the front everywhere is small. In particular, if the entire front corresponds to a constant time function, then 
the algorithm can continue exactly as if this constant-time front were the initial input. To make the entire front or some 
subset of the front coplanar, the algorithm may have to limit the amount of progress at each step so that the vertex 
being pitched is coplanar with some or all of its neighbors. This coplanarity constraint, in addition to causality and 
progress constraints, further reduces the temporal aspect ratio of resulting tetrahedra and may create near-degenerate 
elements. Therefore, it is important to anticipate coplanarity constraints by limiting the amount of progress in the 
current step and thus ensure that both current and future tentpole heights can still be provably bounded from below. 

In this chapter, we show how to make tent pitching less greedy about maximizing the progress at each step so that 
a subset of the front can be made coplanar in the next step. The central idea of the algorithm in this chapter is to trade 
off the amount of progress, i.e., tentpole height, achieved in the current step and that achieved in the next step. This 
algorithm can be used in arbitrary dimensions to make the entire front conform to a target time T > 0. The algorithm 
can also be used to coarsen the front and to retriangulate it by edge contraction or edge dilation by first making a 
subset of the front conform to a local target time. The main feature of the algorithm in this chapter is that it retains the 
progress guarantee up to a constant factor proved for the algorithms in earlier chapters. Thus, we are able to achieve 
target time conformity without significantly sacrificing the worst-case temporal aspect ratio of spacetime elements. 
We believe that the central idea of the algorithm in this chapter will be useful for other purposes such as when the 
spacetime mesh is required to conform to a singular spacetime surface, not necessarily a constant-time plane. 

6.1 Conforming to a target time 

We say that a spacetime mesh Q. conforms to a target time T (or in general to any front 4>) if the portion M x T of the 
constant time plane t ~T is equal to the union of causal facets of Q.. 

In the context of pitching tents, conforming to a target time T means that the height of a tentpole is further 
constrained so that the top of the tentpole does not exceed time T. This additional constraint creates a problem when a 
vertex P = {p, t(/?)) has almost achieved the target time but not quite, i.e., when T{p) < T and T — T{p) is very small. 
In this case, advancing p from T{p) to T will create a tentpole of height T — T{p) which is unacceptably small and 
results in degenerate or poor quality elements. In this section, we show how to incorporate additional constraints into 
our algorithm so that the advancing front conforms to the specified target time T while still proving a lower bound on 
the height of each tentpole. Thus, we retain the quality guarantee of the algorithm within a constant factor. 
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Let 7 e (0, J ] be an additional parameter to our advancing front algorithm. 

The key idea we use is that whenever the front attempts to get too close to the target time T during a single step 
of pitching a local minimum p, we reduce the height of the current tentpole so that in the next step the vertex p will 
achieve the target time T . We call the current step the "penultimate" tent pitching step because the next step is the last 
step in which vertex p is pitched. The tradeoff we need to make is between the height of the tentpoles in the current 
(penultimate) step and in the next (final) step. We do it in a way that in both cases the tentpole height is at least 7 times 
the worst-case guarantee on the height of any tentpole at p. 

Let T denote the front at the current iteration. Let P be a local minimum of T such that z{p) < T. Such a local 
minimum vertex p must exist unless the entire front T has achieved the target time T; in particular, the global minimum 
vertex is a candidate vertex. 

Denote by Wp the minimum distance of p from aff A for every simplex A in (the spatial projection of) the front T 
incident on p. 

Wp:= min {dist(p,aff A) } 

AGlink(p) 

Let (7 = Oimn denote the reciprocal of the global fastest wavespeed. For any vertex p, let hp denote a lower bound 
on the height of the tentpole at p whenever p is pitched. Note that we require hp to be small enough to be independent 
of the front during the step in which p is pitched, i.e., hp can be a function of only the spatial projection of any front. 
Depending on which progress constraint the algorithm uses, i.e., depending on the specific meshing algorithm, hp may 
be a different lower bound; thus, for instance, the algorithm of Chapterl2]guarantees a worst-case tentpole height of 

hp := min{£, 1 — e}wpa. 

The algorithm at each step enforces the following invariant on the front T at the beginning of step /: for every 
vertex /? of T we have: 

either T(p) = T or t(/?) < T - 7/ip (6.1) 

We assume, of course, that the initial front Tq satisfies this invariant, i.e., either the space mesh M is refined enough or 



the target time T is large enough to satisfy Equation 6. 1 



Suppose the algorithm is currently pitching such a local minimum vertex p of T such that T{p) < T. Such a 
candidate vertex must exist because a global minimum of the front is always a candidate unless the entire front is at 
time T. 

Let t'{p) denote the time value of the top of the tentpole at p, maximized to satisfy all causality and progress 
constraints. The new algorithm is allowed to lower the top by choosing a different time value T"{p) for the top of the 
tentpole such that x{p) < x"{p) < T'{p). 

Let 7 be a fixed constant in the range < 7 < 5, a parameter to the algorithm. For the reasons below, 7=5 may 
be a good choice in practice. 

Apply the following rule at each tent pitching step. 

1. Case 1: If i'{p) > T, then choose t"{p) := T as the new top of the tentpole at p. 

2. Case 2: Otherwise, if t'{p) >T — yhp, then choose T"{p) := T — (1 — y)hp as the new top of the tentpole at p. 

3. Case 3: Otherwise, we must have 'l'{p) <T — 7/1^; then choose T"{p) := ^'(p) as the top of the tentpole at p. 
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We prove that the height of the tentpole in the current step is possibly reduced to yhp instead of hp; however, the 
height of the tentpole the next time p is pitched is guaranteed to be at least (1 — Y)hp. The choice of parameter 7 
represents a tradeoff between the current step and the next step. 

As a corollary, the entire front will achieve the specified target time T in a finite number of steps. 

It is easy to see that the algorithm must eventually advance every vertex that has not yet achieved the target time 
value T. We prove next that whenever a vertex p is advanced, then the height of the resulting tentpole at p is at 
least yhp > 0. 



First, we prove by induction on the number of tent pitching steps that the invariant of Equation 6.1 is always 



maintained. Let T, denote the front at the beginning of iteration /. We assume that the initial front To automatically 



satisfies the invariant. Whenever case 1 applies at step /, we have T,+i(/7) — T; therefore, the first part of Equation 6.1 
is satisfied by the new front. Otherwise, if case 2 applies, then T,+i (p) :—T — {l— y)hp <T — yhp because 7 < 1 — 7. 
Finally, if case 3 applies, then we have T,+i (/?) :— i'{p) <T yhp, so the second part of Equation 6. 1 is satisfied by 
the new front. 

Theorem 6.1. Let T be a progressive front. Let y be a parameter in the range (0, j]- Let T be a target time such that 
T > t(/?) + yhp for every vertex p. Then, whenever our new algorithm advances a local minimum vertex p of X, it 
constructs a tentpole at p with height at least yhp. 

Proof We prove that T{p) < T"{p) < t'{p) and that t"{p) — T{p) > yhp. 

Note that t"(/?) = r - (1 - y)hp <T- yhp < T'{p) because 7 < i. 

We show that lowering the top of the tentpole at p from T'{p) to t"(/?) does not create degenerate elements in the 
current step. 

There are three mutually exclusive possibilities that cover all eventual outcomes. 

(a) Suppose case 1 applies in the current step, i.e., T"{p) = T. Since T{p) < T, it must be the case that either case 2 
or case 3 must have applied in the previous step. If case 2 applied in the previous step, then T{p) — T — {l — y)hp, 
so the height of the current tentpole is (1 — y)hp > yhp because 7 < 5. 

This is the only (sub)case for which the claim is tight. 

On the other hand, if case 3 applied in the previous step, then T{p) <T — yhp, so the current tentpole height is 

T~T{p)>yhp. 

(b) Suppose case 2 applies in the current step, i.e., T"{p) = T — (1 — y)hp. Then, the current tentpole height is 

x"{p) - x{p) and we have x"{p) - x{p) ^{T-{\- y)hp) - x{p) = [T - x{p))- (1 - y)hp. 

But from the tent pitcher progress guarantee, we know that t'{p) > T{p) + hp. 

Hence, r - t(/7) > {T -x'{p))+hp. 

We also know that T'{p) < T because case 2 applied instead of case 1. Hence, T — T{p) > hp. 

Therefore, t"(/?) - t(/?) > hp - (1 - y)hp = yhp. 

(c) Suppose case 3 applies in the current step, i.e., t"(/?) = t'{p). Then, from the tent pitcher progress guarantee, 
we know that T'{p) ~ T"{p) > T{p)+hp. 

n 
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As an interesting corollary, we obtain the following result for the linear non-adaptive case. We can construct a 
spacetime mesh with the same progress guarantee that conforms to any progressive front <i>. We can do this by running 
our algorithm backwards in time, starting from <I> as the "initial" front and choosing to conform to the constant time 
plane f = 0. (Conceptually, we can imagine this as meshing the spacetime volume M x [— <I>,0].) It is possible to run 
the algorithm backwards in time with the same progress guarantee because the gradient constraints on the front at each 
step are unchanged if all time values are replaced by their negatives. 

Further, we can mesh the spacetime volume between any two progressive fronts <I>i and <I>2 over the same space 
domain M, as long as there exists a constant time plane t — T between 4>i and <I>2 that has a margin of at least 
yTmin from either of the two fronts. We do this by running the algorithm in two phases: once we run the algorithm 
starting from the front <I>i and conforming to time t = T; then, we run the algorithm starting from the front <I>2 and 
conforming to time t~T. By adjoining the resulting two spacetime meshes, we obtain a spacetime mesh of the volume 
between Oi and <p2- The final spacetime mesh is a weak simplicial complex if the spatial projections of the given two 
fronts <I>i and <I>2 are conforming, i.e., each triangle in one triangulation is the union of one or more triangles of the 
other triangulation. 

Anticipating changes due to adaptive refinement and coarsening For a vertex p, the local spatial 
geometry of the front can change due to adaptive refinement and coarsening or due to some other retriangulation of 
the front; thus, Wp can change and therefore a fixed hp no longer suffices as an adequate lower bound on the tentpole 



height at p. However, we can interpret the invariant of Equation 6. 1 as a constraint (an upper bound) on the spatial size 
of the triangulation in the local neighborhood, star{p), of p, rather than as a lower bound on T. If the initial front Tq 
does not satisfy the invariant, then we refine the front until it does. If coarsening a pair of triangles would violate 
the invariant, then this coarsening is not permitted. Any other retriangulation that would violate the invariant is also 
not permitted. However, since refining the front is a requirement for numerical accuracy, we anticipate all possible 
changes in shape of the front due to refinement in advance and alter the definition of Wp appropriately. This is easy 
because newest vertex bisection generates only a small number of homothetic triangle shapes on the front. 

6.2 Conforming to a target time: a simpler heuristic 

The following is another algorithm for making a vertex p conform to its target time Tp while still guaranteeing a lower 
bound on the height of each tentpole height. The following logic is easier to implement than the previous one, because 
it does not involve the progress guarantee computation and does not involve the additional parameter 7 (even though 
it can be modified to depend on such a parameter). 

Let p be a local minimum of the front t; let T{p) denote its time coordinate before lifting and let T'{p) denote 
the time coordinate of the proposed tentpole at p, maximized subject to causality and progress constraints only. Thus, 
h := t'(/:i) — T{p) denotes the proposed height of the tentpole at p before imposing any target time constraint. 

Let Tp denote the target time for p. We will compute a new time value T"(p) for the tentpole top, instead of T'(p). 

Apply the following rules at each tent pitching step: 



I. Case 1: t'(p) > Tp. Choose t"(/?) :== T, 



p- 



2. Case 2: T'(p)<Tp. If T- t'^;-) </;, then Tp < 2t'(/7) - t(/:>); in this case, choose t"(/7) := (rp + f(/7))/2. 
Otherwise, we have T'{p) <T — h,m this case, let T"{p) :— T'{p) remain. 



First, we prove that lowering the tentpole top is acceptable in the current step. We have f"(/?) = {T — t{p))/2; at 
the same time, we are given <T — t'{p) < h which means 

T-t{p)^T-t'{p)+t'{p)-t{p)^T-t'{p)+h>h 

Hence, we conclude t"{p) = {T - t{p))/2 > h/2. 

Finally, we prove that lowering the tentpole top is acceptable in future steps. We have T '-t"{p) — T—{T — t{p))/2; 
we also know that T — t'{p) < h, so T — t{p) — T — t' (p) + 1' (p) —t{p) — T — t'{p)+h < 2h. Hence, we conclude that 
T-t"{p) <T-2h/2=T-h. 

6.3 Smoothing out large variations in tentpole heights 

We would like to smooth out large variations in the heights of the tentpoles erected at any given vertex p in any two 
consecutive steps. 

Let p he a local minimum of the front t; let T{p) denote its time coordinate before lifting and let T'{p) denote 
the time coordinate of the proposed tentpole at p, maximized subject to causality and progress constraints only. Thus, 
h :— t'{p) — t{p) denotes the proposed height of the tentpole at p before imposing any target time constraint. 

Let 7 > be a parameter. Let hp denote a lower bound on the height of the tentpole at p; thus, 

hp :— min{e, 1 — e}wpa. 

If h> (1 + y)hp, then choose h' ~ {h + hp)/2 as the new tentpole height; otherwise, let h' — h remain the chosen 
tentpole height. Clearly, this choice still retains the progress guarantee of hp. 

In the context of a target time constraint, we apply the above rule only in the last case when the original tentpole 
height is accepted by the logic that aims to achieve the target time Tp. If the tentpole height has already been curtailed 
by a target time constraint, then we do not reduce it any further. 

6.4 Coarsening 

Coarsening is always stated as a request while refinement is a necessity. However, for the purpose of efficiency, it is 
important to coarsen aggressively any parts of the front that are marked as coarsenable so that the number of spacetime 
elements, and thus the total computation time, is reduced. 

Coarsening the front by merging adjacent triangles requires the triangles to be coplanar in spacetime. In the 
previous sections, we gave algorithms to make the entire front conform to a global target time T . We use a similar 
approach to give a new reliable coarsening algorithm. Our new algorithm guarantees that any part of the front that is 
marked as coarsenable will eventually be coarsened. This requires that coarsenable pairs of triangles be made coplanar 
after a finite number of tent pitching steps. Often, this requires that the actual tentpole height at a vertex of the two 
triangles being merged be reduced to achieve the coplanarity constraint. Our new algorithm guarantees that, even with 
the additional coplanarity constraints, the worst-case tentpole height, and hence the minimum temporal aspect ratio 
of elements, is bounded from below. In fact, the minimum tentpole height at a given vertex is guaranteed to be at 
least 7 times the minimum height in the absence of any coplanarity constraints, where y e (0, \] is a parameter to the 
coarsening algorithm. 



Suppose 5 is a subset of triangles that are scheduled for coarsening. Choose a synchronization time T and assign it 
to each vertex incident on a triangle of S. Now, use the algorithm from the previous section to ensure that each of these 
vertices achieves their target time T . When all of them are coplanar with the constant time plane t — T, then perform 
the scheduled coarsening and then remove the scheduling constraint T from all vertices. Recompute a new target time 
for all vertices of the coarser triangles if the coarser triangles are also coarsenable. In fact, we also recompute the 
target time for every vertex of a triangle that gets refined. 

What if multiple coarsenings are scheduled at the same time? In that case, it is possible that some vertex p has 
achieved its scheduled time T2 and is waiting for another vertex i of a sibling triangle to also achieve time 72. Note that 
s need not be a neighbor of p, but s must be a neighbor of an immediate neighbor of p. In such a situation, p cannot 
be pitched even though it is a local minimum because p is waiting for s. However, this is not a problem because we 
can prove that s must eventually achieve time T2. 

If s has no scheduling constraint, then it will eventually be pitched because it will eventually become a local 
minimum vertex. Since s is lower than time T2, it will be pitched to time T2. On the other hand, if s has been scheduled 
to achieve a target time Ti because some other triangle incident on s is being coarsened, then it must be that Ti < 72. 
Hence, by induction, eventually s will achieve time Ti at which step the algorithm is free to pitch s further until it 
eventually achieves time 72. 

Formally, the new algorithm with coarsening constraints is given next. 

6.4.1 New algorithm 

Let (7 = (7min denote the reciprocal of the global fastest wavespeed. 

Given a front T, and a vertex p of the spatial projection of T, denote by w' the minimum distance of p from line qr 
for every triangle pqr incident on p. 

w':= min {dist{p,affqr)} 

Apqr3p 

Let e denote min{e, 1 — e}. Since < e < 1, we have < e < 5. Let h'p denote ew'pG. Note that w^ and hpi depend 
on the current front T,. 

Call a vertex v a leaf vertex if v is an interior vertex incident on four triangles, or v is a boundary vertex incident on 
two triangles, and if every triangle incident on v is a leaf in the refinement hierarchy. Call a leaf vertex v coarsenable 
if level(v) > and if all triangles incident on v are marked coarsenable by the numerical analysis. 

Initialize 

Fix the value of parameter 7 such that < 7 < 5, e.g., J = \- Initialize the front Tq. Every vertex of To is assigned 
level zero. Every triangle of To is assigned a level number of zero. 

Every vertex /? of T (not just a local minimum) is assigned a target time Tp. We maintain the invariant that for 
every front T, for every vertex /? of T we have T{p) < Tp. 

If the problem requires a finite global target time, then let T denote this target time, otherwise let T <— 00. Initially, 
every vertex p of T (not just a local minimum) is assigned a target time Tp ^- T. 

Refine the current front T as many times as needed until the following inequality is satisfied for each vertex p: 

Tp>T{p) + yhp 
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Note that hp is halved after every two levels of newest vertex bisection; hence, for a fixed Tp the refinement process 
will terminate. The refinement process does not change T{p), hence the above inequality gives an upper bound on hp. 

Tent pitching 

At the current iteration of the algorithm, let T,- denote the current front at the beginning of iteration /. 

Step 1: Assign a target time to each vertex (not just local minima). The general idea is that if 
a set of triangles of T is marked as coarsenable by the numerical analysis, then we bring all the five vertices involved 
to a common target time value; then, the triangles are coplanar and the front can be coarsened immediately. 
There are three steps to computing a target time Tp for each vertex p. 

(1) For each vertex p, there is a lower bound Ip on the target time Tp. (The target time Tp is computed in step 3.) A 
lower bound is necessitated in order to maintain the invariant of Equation 6. 1 This lower bound Ip depends on all the 
triangles on the current front incident on p. 

The lower bound Ip on the target time Tp is computed as follows: Ip — Ti{p) + jh'. 

(2) We can coarsen a set of four triangles at once (as long as they are marked for coarsening by the numerical 
analysis). Call this a coarsenable cluster C. Compute a target time Tc for the cluster as the maximum of Ip for each 
vertex p in the cluster. (There will be exactly five vertices incident on the four triangles.) 

If s denotes the degree-4 vertex of the cluster C that will be removed as a result of coarsening cluster C, then the 
vertex s participates in only one coarsenable cluster, the cluster C. Therefore, the target time Tc is also the target 
time Ts assigned to vertex s in step 3. 

(3) For each vertex v, compute the target time Ty as the minimum of the cluster times Tc for each coarsenable 
cluster C to which the vertex v belongs. Since Tc > U- for each cluster C, we conclude that Ty > /„ and, hence, the 



invariant of Equation 6.1 is satisfied. 

Note that a cluster of triangles is coarsenable only if each of the four triangles in the cluster is a leaf, i.e., has not 
been subdivided any further If v does not belong to any coarsenable cluster, then the target time T,, is the global target 
time T (which may be +°o). Note that step (1) looks at all triangles incident on p in order to compute Ip, but step (3) 
looks at only coarsenable triangles incident on p in order to compute Tp. 

After a coarsening takes place, a new lower bound Ip is computed for each vertex p of the two coarser triangles 
that violates the invariant, and a new target time is computed for each coarsenable cluster to which p belongs. 



Let s be a coarsenable vertex (see Figure 6. 1 1 and let C denote the set of four coarsenable triangles incident on s. 



(See figure.) Let V be the set of five vertices of the triangles in S. (In Figure 6. 1 the set V is {p,q, r, s, u}.) 
For each vertex v G V, compute Z,, such that 

iy = Tiiv) + rh[. 

Next, assign a target time T, to the coarsenable vertex s (and thus to C) such that T, := max,.GV /,.. 
Mark each vertex p such that T{p) < Tp as not ready. 

Step 2: Pitch a tent at a local minimum. Let p be an arbitrary local minimum of T such that t(/?) < Tp. In 
other words, select as a candidate for tent pitching only those local minima that have not yet achieved their respective 
target times. At least one candidate vertex must exist; in particular, a global minimum vertex must be a candidate 
vertex unless the entire front has achieved the global target time T. 
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Figure 6.1: Coarsening a set of four triangles in one step. 

Let t'{p) denote the time value of the top of the tentpole at p, maximized to satisfy all causality and progress 
constraints. The new algorithm chooses a possibly different time value T"{p) for the top of the tentpole such that 

T{p) < T"{p) < T'ip). 

Apply the following rule at each tent pitching step. 

1 . Case 1: If t'(p) > Tp, then choose T"(p) := Tp as the new top of the tentpole at p. After pitching p, mark p as 

ready. 

After pitching p, some cluster C of triangles participating in a single coarsening step may become coplanar with 
the corresponding target time Tc- For each coarsenable vertex q neighboring p, if all neighbors of q are ready, 
then all triangles incident on q are coarsened immediately. Note that q is deleted as a result of this coarsening 
step. For each neighbor of q, compute a new target time as described in Step 1 . 

2. Case 2: Otherwise, if t'(/:i) > Tp — yhp, then choose t"(/:)) := Tp — {1 — y)hp as the new top of the tentpole at p. 

3. Case 3: Otherwise, we must have t'(p) < Tp — yhp\ then choose t"{p) :~ T'{p) as the top of the tentpole at p. 

If a patch is rejected, causing refinement of the front, then each new triangle is marked not coarsenable. Whenever 
an edge is bisected as a result of newest vertex bisection, then assign the new vertex (the interior point of the bisected 
edge) a target time equal to the global target time T. 



A cycle in the propagation pattern 

If the pattern of marked vertices is such that it causes a cycle in the propagation pattern (see Figure [6^ , then the only 
reliable way to coarsen any pair of triangles in the set is to coarsen them all simultaneously. Fortunately, if the largest 
edge of each triangle is chosen as the base, then a cyclic marking is very unlikely to occur because it would require 
a cycle of isosceles triangles. If ties between two edges of equal length are broken consistently, for instance using 
lexicographic order, then such a cycle will never occur 



6.4.2 Proof of correctness 

It remains only to prove that the minimum height of each tentpole is bounded. In the previous section, we showed that 
for every vertex p of the front T, the minimum height of every tentpole was at least 7/1' Since there is a minimum 
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Figure 6.2: Coarsening a set of several triangles in one step. 

element size that is always numerically acceptable, there is a global minimum hp that implies a corresponding positive 
lower bound on the minimum height of each element. 

Chapter summary 

Guaranteed coarsening has been an open problem ever since our adaptive meshing result was published ( |ACE+04| |. 
We need to strike the right balance between efficiency and accuracy — the former achieved by adaptive refinement and 
aggressive coarsening, the latter achieved by theoretical guarantees on worst-case element quality. The coarsening al- 
gorithm in this chapter retains the lower bounds on element temporal aspect ratios proved in the absence of coarsening. 



and this represents a big step forward. However, the heuristic of Section 6.2 seems to perform a little better in practice. 



We can imagine a few other heuristics that would attempt to improve element quality and solution accuracy in practice, 
but it seems hard to prove that they will be efficient in terms of guaranteed coarsening of coarsenable regions. 

A new coarsening operation proposed by Jeff Erickson is to delete an interior vertex U of degree-4 with neighbors 
P, Q, R, and S by constructing a patch consisting of two tetrahedra: either the two tetrahedra PQRU and PRSU, or the 
two tetrahedra PQSU and QRSU; the first choice creates two new triangles PQR and PRS on the new front while the 
second choice creates the triangles PQS and QRS. The coarsening operation is performed only when the quality of each 
tetrahedron in the patch exceeds an acceptable threshold. The same idea is used to coarsen the boundary by deleting a 
degree-3 boundary vertex and creating a single-tetrahedron patch. The algorithm extends naturally to deleting a vertex 
u of arbitrary degree; the triangulation after deleting u is obtained, for instance, by contracting the edge up for some 
neighbor p of u. We could prefer contracting the shortest contractible edge incident on u. The alternative coarsening 
algorithm is to pitch possibly inclined tentpoles at u and/or any of its four neighbors, to make triangles PQU and UQR 
coplanar and the triangles PUS and SUR coplanar before merging them pairwise. The former coarsening method does 
not guarantee non-degeneracy of the spacetime elements but uses fewer tetrahedra (two instead of four) and does not 
require that pairs of triangles be coplanar. 
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Chapter 7 

Extensions 



In this chapter, we describe extensions to our algorithms, in the form of newer advancing front operations. We show 
how to incorporate operations that simultaneously advance the front in time as well as smooth the spatial projection 
of the front to improve the spatial aspect ratio of front simplices. Front smoothing will likely increase the amount of 
progress in subsequent tent pitching steps and potentially improve the quality of future spacetime elements. These 
more general operations are potentially useful for tracking moving boundaries. We give a framework that employs the 
new advancing front operations to make local modifications to the front in response to changes in the geometry of the 
domain. We describe our recent work and initial progress on tracking moving domain boundaries. 

Next, we discuss details of implementing our algorithms, such as the choice of various parameters. Our advancing 
front algorithms can be tuned to perform much better in practice than the theoretical performance guarantee by a 
suitable choice of parameters and heuristics. 

7.1 Extensions and heuristics 

In this section, we describe extensions to our advancing front algorithm and propose heuristics to improve the perfor- 
mance in practice. 

7.1.1 Asymmetric cone constraints 

In the presence of anisotropy, waves travel with different speeds in different directions. Therefore, a cone of influence 
whose slope in spatial direction n is equal to the wavespeed in that direction must necessarily be non-circular. The 
resulting cone constraint is asymmetric because it enforces a different gradient constraint in different spatial directions 
in order to satisfy causality. 

Our algorithm can be easily extended to handle asymmetric cone constraints by enforcing the causality constraint 
separately for each spatial direction. Thus, we enforce a constraint on the gradient of every front in every spatial 
direction n to satisfy the causality constraint in the direction n. This will be our approach in higher dimensions as well. 
In the presence of asymmetric cones of influence, we assume implicitly that each slope a is quantified universally for 
every spatial direction n. Causality and progress constraints are interpreted as gradient constraints in each spatial 
direction n. 

We can re-interpret each theorem in the anisotropic setting. We understand every statement to be implicitly quan- 
tified over every spatial direction n. 

For instance, we can restate our theorems for 2D x Time in the anisotropic setting as follows. Let APQR be a 
triangle of the front T with P as one of its lowest vertices. Without loss of generality, assume x{p) < x{q) < T(r). 
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Let e be any real number in the range < e < 1 . Let Gynm denote the global minimum causal slope in any spatial 
direction. Let C7n( A) denote the minimum causal slope in spatial direction n over every point of A, i.e., 

C7„(A):-min{(7„(P)} 

PeA 

1. Let AP'QR denote the triangle of the front t', obtained by pitching the local minimum P, for arbitrary T'{p) in 
the range t(^) < t'(;?) < t(^) + e(7rnindist(/?,aff^r). 

For every spatial direction n e R^, if || V t|^,J < (1 - e)^gr(yn{,P'QR), then l\P'QR satisfies j] Vt']] < Oa^P'QR). 

2. Let AP'QR denote the triangle of the front t', obtained by pitching the local minimum P, for arbitrary T'(p) in 
the range t(^) < T'{p) < T{q) + (1 - e)c7mmmin{|/:>^|, |/7r|}. 

For every spatial direction n e M", AP'QR satisfies 

II V ^\J < {l-e)c^p,an{P'QR) and ||V t'\J < {l-e)c^,„-CTn{P'QR) 

Our algorithm is modified to enforce gradient constraints simultaneously in every spatial direction. In other words, 
we require that for every triangle PQR with local minimum vertex P, the slope of the edge QR is less than or equal to 
(1 — e)^qr<^n{PQR). If the wavespeed is increasing in the future, then the lookahead process will enforce a stronger 
gradient constraint parameterized by the future wavespeed a^. All previous theorems assure positive progress at each 
step when the appropriate gradient constraint is enforced for each spatial direction. The correctness of our algorithm 



requires the anisotropic no-focusing axiom. Axiom L3 which is also sufficient. 

The asymmetric nature of cones of influence complicates the data structures, specifically the bounding cone hi- 
erarchy. While constructing the hierarchy, at each step we need to compute a tight, possibly asymmetric, cone that 
contains two given cones of influence. We have several choice for representing asymmetric cones, such as polyhedral 
cones (the convex hull of vectors in spacetime with the same origin) or cones with elliptical cross-sections. For the 
sake of efficiency, we require each cone to have constant complexity. 

However, a noncircular cone is still a ruled surface. Therefore, both our earlier algorithms, exact and approximate, 
to maximize the progress at each step subject to causality continue to apply without any changes. 

Anisotropic response means that waves propagate faster in one direction than another This can lead to nonlocal 
cone constraints. However, our algorithm already handles nonlocal constraints. 

7.1.2 Non-manifold domains 

So far, we have assumed implicitly that the space domain at time f is a li-dimensional manifold. However, our 
algorithms and all discussions in this thesis extend to non-manifold domains. An example of a non-manifold domain 
in practice would be three thin metal sheets intersecting in a T-junction, or a mechanical assembly consisting of parts 
with complicated intersections. 

Let Mt denote the space domain at time t. The space domain Mt is a subset of Euclidean space of appropriate 
dimensionality such thatM; has a ^/-dimensional simplicial complex. In this case, Mt could be a collection of triangles, 
edges, and vertices forming a 2-dimensional simplicial complex embedded in E^. However, M, need not be a surface; 
for instance, M, could consist of three triangles sharing a common edge. In this case, we say that the space domain has 
dimension d — 2. We measure distances between points in M, by taking the corresponding distances in the ambient 
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Euclidean space which by the nature of the embedding is a lower bound on the geodesic distance between the same 
points. 

7.1.3 Front smoothing 

The progress constraint at the /th step also depends on the geometry of the front in step ; + 1 . As in the previous 
section, our aim is also to ensure that the progress guarantee at the /th step is proportional to local geometry of the 
spatial projection of the front T,. The discussion in this section applies to arbitrary dimensions d>2. 
The amount of progress made by a local minimum p of T, is proportional to 

Wp:= min dist(/7,aff A). 

A€lmk{p) 

If we can increase Wp, then we would get taller tentpoles and therefore fewer spacetime elements overall. In this 
section, we give an algorithm that allows more general operations such as pitching inclined tentpoles and performing 
edge flips as long as causality and all progress constraints are satisfied. Each such operation is desirable if it improves 
the front locally in the sense of increasing the minimum Wp for any vertex p in the local neighborhood. 

Inclined tentpoles 

Mesh smoothing in the neighborhood of a vertex p is achieved by pitching an inclined tentpole at p from P= {p, T{p)) 
to P" = {p' ,x'{p')). Advancing p in time as well as moving it in space can be thought of as a combination of two 
successive operations: (i) move p in space to p' , and then (ii) advance p' in time from P' ~ {p',T{p')) to P" = 
(/?', t' (/?')). Each of the spacetime elements in the resulting patch shares the tentpole PP" and has both a causal inflow 
facet on the old front T as well as a causal outflow facet on the new front t'. 

We assume that the direction of the incUned tentpole PP" is fixed, for instance in an earlier iteration of the algo- 
rithm. To perform Laplacian smoothing, the direction of PP" could be chosen to move p closer to the centroid of its 
neighbors in the spatial projection of T,. 

Let f denote the intermediate front obtained from T after the first step. If f is progressive and if P' is a local 
minimum of f, then the final front T,+i = advance(f,/?',Af) after pitching P' = {p' ,f{p')) to P" = (/?',T,+i (/?')) is 
progressive, for some At > 0. It follows by induction on the number of tent pitching steps that every front constructed 
by the algorithm is valid. 

The length of the inclined tentpole PP" is constrained so that every triangle P"Q' fr is causal, i.e., not too steep. 
See Figure [7T| 

Anticipating changing geometry 

The primary challenge we face is how to modify the progress constraints in response to changing geometry of the 



spatial projection of the front. The need to strengthen the progress constraint is clear from Lemma 2.4 — if ^{pqr) 



decreases because the obtuse angle of l\pqr becomes more obtuse (Figure 7.2 1, then positive progress cannot be 
guaranteed. The factor {pqr) that appears in the progress constraint of Definition 2.6 must be changed to anticipate 
any potential increase in the largest obtuse angle of any triangle of the front. 

We enforce a lower bound on the smallest acute angle and on the largest obtuse angle in the spatial projection of 



any front. Fix (p^{^ > such that 0nun is less than <p{pqr) for any triangle pqr of the front T. (Recall Definition 2.5 ) 
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Figure 7.1: The length of an inclined tentpole is limited so that new front triangles are causal. 




Figure 7.2: Moving a vertex can introduce or increase an obtuse angle on the front. 



We will enforce that sin 6 > (j)jnm for every angle of every triangle on the front constructed at each step. We are now 
able to strengthen the progress constraint to anticipate the changing geometry of the front as long as the front does not 
violate the following angle bound. 

Definition 7.1 (Angle bound (pmm)- W^ s<^y that a triangle Aabc satisfies the angle bound (^min if the minimum sine of 
any angle of Aabc is at least (pmin- 

In many ways, the choice of (/)n,in is like the estimate of the global maximum wavespeed. We will exert a lot of 
effort in Chapter l3] to obtain better estimates of future wavespeed than the global maximum. However, for arbitrary 
changes to the front due to spatial motion, we do not know yet how to adaptively obtain a better estimate of how 
(j) (pqr) changes. For now, we have to be content with a global bound on the smallest acute angle as well as the largest 
obtuse angle allowed in the spatial projection of any front. Techniques similar to those in Chapter[3]may prove useful 
to adaptively estimate the changing geometry due to motion. 

Currently, we do not know how to exploit the fact that the algorithm may freely choose the tentpole direction. 

We strengthen the progress constraint by replacing, in Definition 2.6 every occurrence of by (j)mm, where ^^in is 



a parameter Since our new progress constraint is stronger. Lemma 2.4 still holds. The price we pay for strengthening 
the progress constraint is a smaller worst-case progress guarantee, smaller by a factor of (pmin, due to the fact that the 
front in the next step has to satisfy a stronger gradient constraint. 
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Definition 7.2 (Progress constraint a). Let PQR be an arbitrary triangle of a front T. We say that the triangle PQR 
satisfies progress constraint a if and only if for every edge, say qr, we have 

II V AqrW — 1— I < (1 -e)CT0„™ 



In arbitrary dimensions, fix 0rnin > such that ^^^ is less than ^{popipi- ■ -Pk) for any A;-dimensional simplex 



PoPiPi- ■ -Pk of the front T. (Recall Definition 2. 11 ) 



Consider an arbitrary ^-dimensional face popip2- ■ -Pk incident on po. 

Definition 7.3 (Progress constraint). We say that a simplex PoA^2- • Pk of the front T satisfies the progress constraint 
(7 if and only if for every facet At where < i < k we have || V t|^.|| < (1 — e) (7(j),„m- 



Note that the progress constraints of Definition 7.3 and Definition 7.3 limit the gradient of every facet of a front 



simplex uniformly, unlike Definition 2.6 and Definition 2.13 which limits the gradient of only highest facets of each 
front simplex. 

Let p he a local minimum of T. Suppose we advance p in time from T{p) to f{p) but no higher than its lowest 
neighbor, i.e., f{p) < T{q) for every neighbor q of p. Additionally, we move p in space to /?' but only so far that 
Ap'qr also satisfies the angle bound 0min- By the monotonicity of the progress constraint, we know that the new 
triangle P'QR where P' = {p',f{p)) satisfies the new progress constraint a as long as APQR also satisfies the new 
progress constraint a. We have thus proved the following lemma. 



Lemma 7.4. If a front X satisfies the progress constraint O of Definition 7.2 and if p is a local minimum of T, then 



the new front f satisfies the progress constraint (J of Definition 2.6 where f is obtained by moving P — (p,t(/))) to 
P' = {p',f(p)) such that P' is also a local minimum of t, p' is in the kernel of star(p), and f satisfies angle bound 

ymin- 

We restrict the motion of p so that for every triangle pqr incident on p, the smallest angle of Ap'qr is bounded from 
below. To build on earlier work which guarantees that p' can make finite positive progress if p' is a local minimum, 
we also need to guarantee that f is progressive and that p' is a local minimum of the front f. It is easy to see that, as 
long as p' stays a local minimum, the front f is progressive if T is progressive; this is because as long as p' is a local 
minimum of f we have ||V f|^|| < ||V t|^|| for every simplex A. Therefore, p' must satisfy x{p') < z{q) for every 
neighbor q of p, where p is a local minimum of T. 



The allowed region for p' such that Ap'qr satisfies the angle bound 0niin is shown in Figure 7.3 when a single 
triangle pqr is considered in isolation where T{p) < T{q) < T(r). When all triangles incident on p are taken together 
and we know that the intersection of their allowed regions is nonempty because it contains at least p. 

Lemma 7.5. If APQR of the front T with x{p) < x{q) < T{r) satisfies ||Vt| || < (1 — e)0„j,„(7, then for every 
At e [0, (1 — e)^„„„dist(/:i,affq'r)(j] the front x' — advance {T,p, At) satisfies || V t'\ \\ < (1 — £)^min<y- 



Proof of Lemma\73\ By Definition 2.6 the front t' satisfies progress constraint a if and only if 



max{ ||V t'|„ JIJV t'I „.||,||V t'I l| } < (1 -e)0„in(T 



\pq 



Since advancing p in time does not change the time function along qr, we have || V t'| .|| = || V t| || < (1 — e)(j)mm<^ 




Figure 7.3: The allowed region for p' is the shaded region in the top: p' is restricted to the allowed region to ensure 
that minjsin Z/?'^r, sin Z/?'r^, sin Zr/?'^} > 0mm. 



Since APQR satisfies progress constraint a, we have 



max{ II V Tl J|, ||V T|^,J, ||V t|, l| } < (1 - e)</.mmCT 



\qr 



rp\ 



Also, /? is a lowest vertex of t; so t(/?) < min{T(^), T(r)}. 

Consider the constraint || V t'| || < (1 — e)0minO'- As long as T'{p) < T{q), we have ||V t'| || < ||V t| || < 
(1 — e)0minO'- Similarly, as long as x'{p) < T(r), the constraint || V t'| || < (1 — e)0minO' is automatically satisfied 
because ||VT'|,.p|| < ||VT|,p|| < (1 -e)(/)mm(T. 

WhenT'(/:>) > t(^), the constraint ||V t'|^,^J| < (1 -e)0minO' is equivalent to t'(/?) < t(^) + (1 -e)0mmO'|/?^|. We 
have 



r'ip) < i{p) + (1 - e)^mmOdp 

<T(^) + (l-e)(/)minCTc/p 
<T(^) + (l-e)(^minC7|Ml 

where the last inequality follows because dp — dist(/7, ?Lffqr) < \pq\. 

Similai'ly, whenT'(/:i) > T(r), the constraint ||V t'| || < (1 —e)0minO' is equivalent to T'{p) < T(r) + (1 — e)0minO'|/"'| 
We have 

t'(p) <T{p) + {l-e)(^ynm(ydp 
<T(r) + (l-e)0minC^^p 
<T(r) + (l-e)(^minCj|/"-| 

where the last inequality follows because dp — dist(/?, aff^r) < \pr\. D 

We have thus shown that as long as the motion of p to p' is constrained to the allowed region at each step, the 
progress made by pitching p is positive and bounded away from zero. 
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Figure 7.4: An edge is flipped by inserting a single tetrahedron. 

Edge flips 

Another front advancing operation is to insert a single tetrahedron PQRS with inflow facets APQR and APRS and 



outflow facets APQS and AQRS; see Figure 7.4 If the volume of the tetrahedron PQRS is positive, this operation 
advances the front. If the solution within PQRS is sufficiently accurate, the single-element patch is accepted. 

Let T and t' denote the fronts before and after this operation. We require that the front t' be causal and progressive 
to guarantee sufficient progress when a tent is pitched in the next step. 

Triangles pqr and prs belong to the spatial projection of T, and triangles pqs and qrs belong to the spatial projection 
of t'. Therefore, the operation results in an edge flip, replacing the diagonal pr of the quadrilateral pqrs with the 
diagonal qs; of course, quadrilateral pqrs must be convex. 

Thus, the edge flip replacing pr by qs is allowed only if pqrs are in convex position in the spatial projection, 
tetrahedron PQRS has positive volume, and all four facets of the tetrahedron are causal and progressive. 

The edge flip is desirable only if the resulting tetrahedron is non-degenerate and has acceptable quality, and if it 
improves or does not worsen the spatial aspect ratio of front triangles, i.e., 

min{(j){pqs),(j){qrs)} > min{(j){pqr),(j){prs)} 

A similar flip operation could be defined in higher dimensions. It is not clear whether the resulting patch, which 
must sometimes consist of more than one spacetime element, satisfies all the objectives such as, for instance, the 
requirement that every element have a causal outflow facet. 
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7.1.4 Heuristics 

In this section, we describe more general meshing operations, in addition to tent pitching, performed as long as each 
operation alters the front only in ways such that the new front is (/;,/) -progressive for the chosen parameters h and I. 
We label these new operations heuristics because we cannot yet provide a guarantee that each such operation will 
be possible when it becomes necessary or even when it is desirable. These enhancements represent a step towards 
an adaptive advancing front algorithm for tracking spacetime features, such as moving boundaries, by aligning mesh 
facets along these features. Developing such a boundary-tracking meshing algorithm is still work in progress. We 
cannot yet guarantee that any of the operations mentioned can always be performed in response to an urgent need to 
improve the quality of the front or to guarantee progress at all. However, we are able to give sufficient and admittedly 
strong gradient constraints under which any front can be retriangulated, subject to the angle bound, so as to achieve 
the generalized front modification operations of this section. 

Note that newest vertex bisection does not improve the quality of front triangles beyond a very limited amount. 
For instance, the largest obtuse angle in a triangle also occurs in its grandchildren. The more general operations on the 
other hand can improve the quality of front triangles, and therefore the worst-case progress guarantee. Our stronger 
progress constraints are sufficient gradient constraints under which a limited amount of improvement of the front can 
be performed whenever permitted by the constraints. 

Our earlier work on adaptivity was restricted to adapting the mesh resolution to an error indicator; now, we are 
concerned also with adapting the mesh quality, in terms of worst-case temporal aspect ratio of spacetime tetrahedra. 

1 . Generalized bisection Our algorithm refines a triangle on the front by bisecting an edge the triangle at a point 
other than its midpoint, as long as the resulting smaller triangles satisfy the so-called angle bound. 

The bisection direction for each triangle must be chosen so that repeated bisection retains the key property that 
each descendant of a triangle PQR has its edges parallel to three of the four directions defined by APQR and its 
apex — the three edges PQ, QR, and PR plus the bisector PS for the appropriately chosen interior point s on the 
base qr. 

Adjacent triangles must be bisected compatibly, i.e., adjacent triangle must agree on how to bisect their common 
edge so that refinement maintains a triangulated front. 

The front is coarsened by undoing previous refinement. 

The new bisection method means, for instance, that the largest angle in any triangle on each front can be bisected. 

To ensure that the front can be advanced successfully after the bisection, each bisection is permitted only if the 
front after bisection is progressive. Note that this condition also limits the choice of direction along which to 
bisect a given triangle. 

2. Mesh smoothing: Advancing a vertex p in time takes P = {p,T) to P' = (/?,t') where t' > T. Pitching an 
inclined tentpole at p takes P = {p, t) to P' = (z?', t') which simultaneously changes the spatial projection of the 
front, taking p to p' such that p' is in the kernel of star(/5). This achieves a limited amount of smoothing of the 
front, also known as r-refinement. 

With the stronger progress constraints, we are able to allow inclined tentpoles subject to the angle bound. As 
long as the slope of each tentpole is subsonic, i.e., greater than the local causal slope, we are still able to prove 
a positive lower bound on the temporal aspect ratio of each spacetime tetrahedron. 
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3. Flipping an edge: We relax the assumption that the spatial projection of each front is a hierarchical refinement 
of the initial space mesh. Flipping an interior edge e ^ {p,q) replaces the diagonal e of a convex quadrilateral 
with the other diagonal e' = {r,s). An edge flip in the spatial projection of a front T corresponds to adding a 
single tetrahedron PRQS to the spacetime mesh; this edge flip can be performed only if the resulting tetrahedron 
has bounded positive volume. 

With the stronger progress constraints, we are able to allow flipping of interior edges as long as all triangles 
of every front satisfy the angle bound. An edge flip is desirable only if the resulting tetrahedron has bounded 
temporal aspect ratio. The algorithm is free to choose the marked vertex of each new front triangle created as a 
result of an edge flip. 

By construction, since each of the above operations is performed only when permitted by the progress constraints, 



the new front t' after the operation satisfies the angle bound and is (/i,/)-progressive. See Lemma 4.8 

Another operation that changes both the spatial projection of a front T as well as the time function is the following. 
Let /? be a local minimum of T. The new progress constraint limits the gradient of each triangle APQR in star(/9) in 
four different directions. Suppose we advance p in time from T{p) to f{p) but no higher than its lowest neighbor, i.e., 
i{p) < T(g') for every neighbor q of p. Additionally, we move p in space to p' but only so far that Ap'qr also satisfies 



the angle bound (j>mm- Then, by Lemma 7.4 we know that the new front f also satisfies the new progress constraint a. 
Therefore, the algorithm is guaranteed to be able to pitch P, a local minimum of the front f, by a positive amount to 
make progress in time to get the new front t' with P' = (/?', t' (/?')). 

7.2 Tracking moving boundaries 

In many applications, the geometry or even the topology of the domain changes over time. For instance, combustion 
of solid rocket fuel eats away at the boundary where the fuel burns. The mesh must adapt by changing the size of 
mesh elements, or by varying the placement of mesh features, or both. How to dynamically re-mesh in response to 
evolving geometry is a well-studied problem jABG^OOt . With other dynamic meshing algorithm, global remeshing 
of the domain is performed whenever the mesh gets too distorted due to motion; at intermediate steps, smoothing or 
other incremental changes are interleaved with the solution. 

In this section, we restrict our attention to meshing in response to changes in the geometry of the domain due 
to boundary motion. We assume that the topology of the domain does not change, which means it is sufficient to 
consider each connected component of the original space mesh separately. In our context, we always remesh locally 
by anticipating and adapting to changes in the shapes of front triangles. An external or an internal boundary (such as 
a material interface) moving with time describes a surface in spacetime. Similarly, shocks describe singular surfaces 
in spacetime. The solution jumps discontinuously across the shock surface. Tracking the boundary means aligning 
vertices, edges, and facets of our tetrahedral mesh with the corresponding surface (approximately). Boundary tracking 
combined with SDG methods implies a highly accurate resolution of discontinuities in the solution field. 

We incorporate spatial motion over time in our advancing front framework by assuming either a spatial velocity 
vector or a direction in spacetime for every vertex of every front. The tentpole direction can be inferred from either 
input. Boundary vertices have the tentpole direction assigned to them by the numerical analysis, which is aware of the 
evolving nature of the domain. For all other vertices, the meshing algorithm is free to choose any subsonic tentpole 
direction. For instance, for the purpose of Laplacian smoothing, we choose for every interior vertex p a direction 
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vector in space pointing towards the centroid of the neighbors of p. The slope of the tentpole at p is chosen to be at 
least a constant factor greater than the reciprocal of the maximum wavespeed anywhere in star(;5). 

We modify the basic framework of our interleaved mesh construction and solution procedure to recompute vertex 
velocities along with or as part of the solution. 

Our motion model is sufficient to capture boundary motion that does not alter the topology of the underlying 
space domain over time — ^boundary vertices remain boundary vertices and interior vertices remain interior vertices 
throughout. This model therefore seems insufficient to capture splitting apart of the domain into more components, 
creation of tunnels, and other ways that increase the topological complexity of the domain such as creation of new 
non-manifold features by folding. 

In previous chapters, we were meshing the spacetime domain, a prefix of M x [0,°°), for the space domain M that 
did not change over time. In this section, our spacetime domain of interest, i2, is the volume swept through time by the 
space domain M, at time t. Since different parts of a front T may have different time coordinates, the spatial projection 
of T need not correspond to the space domain at any time t. 

7.2.1 Tracking boundaries in ID x Time 

Meshing in ID x Time presents fewer challenges because the front at each step is constrained only by causality. We 
are free to modify the front at each step, such as by refinement or by smoothing, as long as the new front is causal. 

Front advancing operations 

Boundary vertices are those that have velocities prescribed by the numerical analysis. We assume that all prescribed 
velocities are subsonic. All other vertices are interior vertices. The velocity of an interior vertex is computed as the 
average of the velocity vectors of its neighbors. 

Our algorithm in ID x Time only performs tent pitching operations, advancing a local minimum vertex at each 
step. The major challenge is prioritizing which vertex to pitch at each step. Since boundary vertices have prescribed 
motion, we have to adapt the geometry of the rest of the front, especially near the boundary, to anticipate the necessary 
change of geometry on the boundary. 

As a first approximation, we choose in our new algorithm to prefer pitching an interior vertex over a boundary 
vertex. 

The goal of our algorithm is to ensure a lower bound on the temporal aspect ratio of spacetime triangles while 
simultaneously satisfying a prescribed global lower bound Zmin on the length of the spatial projection of all spacetime 
triangles. 

Let ;? be a local minimum vertex chosen by the algorithm to pitch next. Let ^ be a neighboring vertex and let r be 
the neighbor of q other than p if it exists. The tentpole at p is inclined and its slope with respect to the space domain 
is equal to the reciprocal of the speed of p. 

The height of the tentpole at p is constrained in three different ways. Let he denote the supremum of the height of 
the tentpole PP' such that P'Q is causal. Let hi denote the height of the tentpole PP' such that P'QR are collinear in 
spacetime. Let h„ denote the height of the tentpole PP' such that the length of the spatial projection of P'Q is equal to 
Imm where /min is the minimum allowed length of the spatial projection of any segment of any front. 

The algorithm chooses the final height of the tentpole depending on the following cases. The algorithm always 
satisfies the causality constraint; hence, the final height h of the tentpole is always less than he- 

Let 5 be a positive quantity, < 5 < 1, large enough to avoid numerical difficulties. Hence, (1 — 5)hc < h^. 
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Note that hj < he always because hi > he contradicts the causaHty of QR. 
First, consider the situation h„ < (1 — 5)hc. 

Case 1 Consider first the case h[,h„ < (1 — 8)hc < he- 

1.1. Case hi < h^ < (1 — S)hc: Choose h :— h„. If the resulting patch is rejected, then bisect pq. 

1.2. Case h„ < hi < {\ — 5)hc: We have no choice but to choose h :~ h„. If the resulting patch is rejected, then 
bisect pq. 

Case 2 Next, consider the case h^ < {I — 5)hc < hi < he- 

2.1. Case /Jw < (1 — d)hc < hi: We have no choice but to choose h ;= h„. If the resulting patch is rejected, then 
choose h :— ^h^ if ^h^ < /!,v, otherwise bisect pq. 

Next, consider the situation h„ > (1 — 5)hc. 



Case 3 Thirdly, consider the case hi < {I — 5)hc < h^. We do not distinguish between the two subcases (i) /i». < he, 
and (ii) h^ < h„. 

3.1. Case hi<{\ — 8)hc < h„: Choose h := he if h^. ~ he is sufficiently large, otherwise choose h := jhe- 

Case 4 Finally, consider the case (1 — 5)he < h„,hi. 

4.1. Case (1 — 5)he <hi<hw: We do not distinguish between the two subcases (i)/z; <hw <he, and (ii) hi <he <h„. 
We have no choice but to choose h:= (1 — 8)he; however, if h„ — he is too small, then choose h :— ^he- 

4.2. Case (1 — 5)he < h„ < hi: We do not distinguish between the two subcases (i) h„ <hi< he, and (ii) hw<he< hi. 
We have no choice but to choose h:= {1 — d)he; however, if h„ — he is too small, then choose h :— jhe- 

Proofs 

We would like to prove that the algorithm makes positive progress at each step while still guaranteeing causality and a 
lower bound of ly^m spatial projection length. In cases 3 and 4, the height of the tentpole is at least a constant fraction 
of he- We claim that in any sequence of tent pitching steps, cases 1 and 2 apply only a finite number of times before 
either case 3 or case 4 applies in the next step. Clearly, causality is always maintained because the final height h always 
satisfies h < {1 — d)he < he- It is easy to see that all four cases together ensure that the spatial projection of every 
segment on each front has length at least Zmin- Hence, all requirements are satisfied. 

Note that cases 1 and 2 deal with the situation /z„, < {l — 5)he. If the resulting patch is rejected, then the segment /5^ 
is bisected which means that the causality constraint /ij, that applies in the next step is smaller than he by a finite amount. 
Specifically, h'e — he = j\en{pq)a. On the other hand, the constraint h^ also decreases to h[^, where h'^ — h„ — 'min<7. 
Therefore, if l^i„ < len{pq)/2, then the constraint he decreases faster than the constraint h„ with repeated applications 
of either case 1 or case 2. Note that len{pq) also halves with each iteration. Therefore, that after a bounded number of 
applications of cases 1 or 2, it must be the case that either case 3 or case 4 applies. 

Thus, we can show that the front advances after a finite number of steps; however, we haven't shown yet that the 
front achieves a target time T in a finite number of steps, unless we also show that the front is refined only a finite 
number of times before being coarsened because h:= hi is chosen at some step. 
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Coarsening 

Let P be an arbitrary local minimum of a causal front T. Let Q and R be neighbors of P. We claim that in the next step 
in which P is pitched to create the new causal front t', it is possible to coarsen and delete at least one of the segments 
PQ or PR or both, without creating degenerate spacetime triangles. 

Without loss of generality, assume T(r) > t(^). First, consider the case that T(r) > T{p) +A for some sufficiently 
large A > 0. Then, in the next step, P can be pitched to P' such that Q — P' — R are collinear in spacetime and 
T'{p) > T{p)+At Where Af > g^ A. 

This suggests that we need to bound the ratio of the spatial lengths of any two adjacent segments on the front. For 
instance, if S^ > |, then bisect pq. Equivalently, if we start with an initial front that satisfies this balance condition 
then the balance can be maintained by propagating refinement of a segment to neighboring segments. 

Next, consider the remaining case that T{q),T{r) ~ T{p)- Then, P can be pitched to P' such that either P'Q or P^R 
is tight against the causality constraint, i.e., has slope equal to (1 — e)cT. This reduces the number of non-strict local 
minima. Specifically, it reduces the number of local minima for which both the incident edges are not tight against 
the causality constraint. By induction on the number of such minima, it follows that eventually there will be a local 
minimum vertex such that at least one of its incident edges has slope equal to (1 — e)(7. (Such a local minimum may 
be a boundary vertex and therefore also a strict local minimum.) In this situation, the first case above will apply and it 
will result in coarsening. 

Guaranteed coarsening In Chapter l6] we give an algorithm that will guarantee coarsening while still main- 
taining lower bounds on the minimum tentpole height. Coarsening a portion of the front takes some time before that 
portion can be made coplanar and coarsenable. As long as the motion of the boundary permits the additional time 
taken to coarsen reliably, for instance if the boundary moves very slowly, we can guarantee that the segments near the 
front will also satisfy the lower bound on the size of their spatial projections. 

7.2.2 Tracking boundaries in 2D x Time 

In this section, we describe an algorithm in 2D x Time that uses the following advancing front operations to track 
moving domain boundaries — tent pitching, smoothing, edge bisection, vertex deletion, and edge flips. Our algorithm 
does not create inverted elements (i.e., tetrahedra with negative volume). We cannot yet guarantee that the algorithm 
will always make progress by creating non-degenerate elements. In the future, we could use a similar algorithm for 
tracking shock fronts. Each operation changes the triangulation of the front. Most operations also advance the front 
and the solution. 

Meshing in 2D x Time to track spacetime features is challenging because of the need to satisfy not only causality 
but also progress constraints at each step. The general framework of our algorithm in 2D x Time is to prefer operations 
that increase the spatial quality of front triangles without introducing refinement whenever possible. We devised 
a policy to assign priorities to various front advancing operations with the aim of proving correctness as done for 
Delaunay refinement algorithms ( IEdeOItlCD03llCP03l l. 

First, we outline the various meshing operations, in addition to the usual tent pitching operation which advances a 
single vertex in time. 

Pitching inclined TENTPOLES We advance P to P' along the given tentpole direction in spacetime, creating 
a tentpole PP'. The tentpole direction is the one computed by the solver; when the boundary is moving the tentpole 
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is inclined with respect to the vertical time axis. The amount of progress, i.e., the length of the inclined tentpole, is 



constrained to prevent inverted elements and is limited by causality. See Figure 7.1 



Smoothing The meshing algorithm can freely choose the tentpole direction for other vertices such as those created 
by refinement in the interior. We choose the tentpole PP' to smooth the front and improve the spatial aspect ratio of 
front triangles. Alternatively, we choose the tentpole PP' to align some implicit facets with a spacetime surface. 

Edge bisection Edge bisection is the same operation as the newest vertex refinement described in Chapter |4] 
Bisecting triangles on the current front decreases the size of future spacetime elements. As long as the vertex incident 
on the largest angle is chosen as the apex of the triangle, newest vertex refinement also decreases the largest angle 
which increases the likely progress in the next step. 

Vertex deletion The front is coarsened by deleting a vertex when permitted by the numerical error estimate. 
Vertex deletion is achieved by inserting a patch of one or more tetrahedra and simultaneously advancing the front. 

Edge flips Flipping an edge of the front triangulation can increase the aspect ratio of front triangles and improve 
the quality of future spacetime tetrahedra. An edge flip of the edge pr, a diagonal of the convex quadrilateral pqrs, is 
achieved by inserting a single tetrahedron PQRS as long as the tetrahedron has positive volume and sufficient quality. 
See Figure [74] 

A framework in 2D x Time 

Choose the applicable operation that is first in the sequence listed below. If a patch created by one of the operations 
is rejected, then proceed to the next operation in the sequence to try to create an acceptable patch. If there are several 
candidates for where to perform a given operation, choose a candidate that makes it most likely that in the next step 
an operation earlier in the sequence will become possible. If there is more than one candidate operation, then perform 
these independent operations simultaneously in parallel whenever possible. 

1 . Coarsen the front by deleting a vertex u 

2. Flip an edge of the front if it improves the spatial aspect ratio 

3. Pitch an interior vertex to smooth local triangulation 

4. Pitch a boundary vertex in prescribed direction 

In general, coarsening operations, or operations that decrease the number of triangles on the front, are preferred 
over those that refine the front because coarsening is harder to achieve in practice and to guarantee in theory. The front 
will still be at least as refined as required by the numerical error indicator. 

Additional heuristics include preferring to pitch local minima in the interior of the domain over those on a moving 
boundary. For an interior vertex, we choose a tentpole direction that is the average of the tentpole directions computed 
for its neighbors. If an obtuse angle gets too big, we flip or bisect the opposite edge. We bisect an edge if its endpoints 
have very different velocities; we choose an average velocity for the new midpoint. 



106 



Each operation is performed only when allowed by the DG solver. We adapt the mesh in response to numerical 
error estimates as before; in addition, we allow the algorithm to refine and coarsen the front in response to changing 
local feature size. 

We speculate that the above policy for choosing the various front advancing operations and prioritizing among 
them is an acceptable first approximation for very simple classes of motion. However, a more sophisticated policy 
needs to be devised. Note that the algorithm incurs the additional complexity of choosing which operation to perform 
at each step. Also, the more complicated the policy, the harder it is to parallelize the meshing algorithm. 

We described a framework using a wider vocabulary of advancing front operations to mesh in spacetime. Even 
without motion, we could use our algorithm to mesh static curved features adaptively. Theorems are lacking, except in 
very restricted cases, such as IDxTime, and uniform translation of the entire domain. It remains a very interesting and 
challenging open problem to devise a provably correct algorithm even for a restricted but still useful class of motion. 

7.3 Practical issues 

There are several choices to be made in practice when implementing our meshing algorithms. 

Choice of parameter e The value of the parameter e in the range < e < 1 controls how greedy the algorithm 
is. Values of £ close to 1 strengthen the progress constraint and therefore potentially limit the height of tentpoles. 
Choosing a very small values of e weakens the progress constraint. It is reasonable to expect that in practice smaller 
values of £ mean that most tentpole heights are limited by causaUty and therefore the average tentpole height is larger. 
However, the theoretical guarantee on tentpole height due to the causality constraint is proportional to e. It is therefore 
important to choose an appropriate value of £, perhaps specific to each problem instance, that achieves both adequate 
minimum tentpole height while maximizing the average tentpole height. A value of e around the middle of its allowed 
range, i.e., £ ~ j, is likely to be a good choice in practice. 

Choice of LOOKAHEAD parameters The horizon parameter h and the adaptive lookahead parameter / (Chap- 
ter l5]l can be fixed a priori or can be increased adaptively if greater progress is desired. Arguably, the algorithm for 
general (/i,Z) is complicated to implement. In practice, depending on the underlying physics, it may suffice to fix the 
parameters at a small value, say h — l— 1 . 

Choice of heuristic Tent Pitcher is not restricted to pitching only local minima, even though the theoretical 
progress guarantee applies only to the minimum progress made by advancing an arbitrary local minimum. To increase 
the degree of parallelism of the algorithm, it is beneficial to pitch any vertex where some positive progress is guaranteed 
even though it might not be the choice that guarantees the most progress at the current step. 

Preliminary experiments suggest that the quality and overall efficiency of the solution strategy depend on the choice 
of which local minimum vertex to advance at each step, whenever more than one candidate is available or whenever 
other operations are also available. All algorithms so far leave this as a free choice to be decided heuristically. We 
strongly recommend making empirical studies of how to exploit different heuristics to generate better quality and more 
efficient meshes. 

Data structures for the front A significant advantage of the advancing front approach is that we need to 
store only the front and a collection of unsolved patches. The front is (i-dimensional, the same as the space domain. The 
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unsolved patches are transient and can be discarded as soon as they are solved and accepted. We recommend separate 
representations of solved elements and unsolved patches because this would allow the solver to be oblivious of how 
the inflow information for a unsolved patch was computed — whether by a previous iteration of the same algorithm, by 
initial and boundary conditions, or by some entirely different solution procedure. Standard data structures JWei85l l are 
available to represent the front in 2D x Time. The adjacency information between facets of the front and facets of the 
spacetime elements stores either a many-one or a one-many relation, because of the weak simplicial complex property, 
which can be represented by pointers or arrays. Thus, our algorithms truly require only very basic data structures 
capable of representing a weakly conforming unstructured simplicial mesh. 

Maintaining a TRIANGULATION When we refine one triangle on a 2D front, we may be forced to refine other 



nearby triangles in order to maintain a conforming triangulation of the front. See Figure 4.3 The propagation path 
touches every triangle in the worst case, but in practice, the propagation path usually has small constant length and 
does not involve loops. This is especially the case if the largest angle in each triangle is used to choose the marked 
vertex. 

The algorithm need not wait for the entire propagation path to be executed before proceeding with the next ad- 
vancing step. In fact, only the initial bisection was required by the numerical error indicator; all further bisections are 
merely artifacts of the algorithm because the algorithm requires the front to be conforming. Therefore, we employ 
a lazy propagation scheme that delays propagating refinement to adjacent triangles unless absolutely necessary. This 
also reduces the need to synchronize and communicate across processors in the parallel setting. 

Lazy propagation splits the adjacent triangle but does not propagate further immediately. Instead, the propagation 
is delayed until the adjacent triangle needs to be refined or pitched. In the interim, the triangulation may consist of 
transient triangles. Eventually, the final result is exactly the same as if the entire propagation path was executed in one 
step. 

Coarsening: A coarsening step that nullifies the last refinement in a refinement propagation path usually does so 



by deleting a degree-4 vertex. An exception to this rule is shown in Figure 4.3 where the last refinement cannot be 
nullified by deleting a degree-4 vertex because all vertices have degree greater than 4. However, a simple edge flip of 
one of the edges involved in the penultimate bisection produces a degree-4 vertex that can be deleted to remove the 
last refinement. Thus, five triangles are involved in the coarsening step instead of just four. Note that this case occurs 
only when a refinement propagation contains a loop. This is not the case when the largest angle of each triangle is 
marked and ties are broken suitably. 

Bounding cone hierarchy Non-constant wavespeed due to nonlinear response means that the strictest cone 
constraint that limits the height of a tentpole can belong to a point of the front arbitrarily far away. It is expensive to 
examine all cone constraints, one for each triangle of the front, to determine the most limiting cone. We adopt standard 
techniques used in computational geometry and collision detection in robot motion planning (ILat91l l. At each step, we 
wish to solve the following optimization problem: maximize the height of the tentpole subject to all cone constraints. 
We group all cone constraints into a hierarchy. Specifically, we use a bounding cone hierarchy to efficiently compute 
a hierarchical approximation of the actual wavespeed at a point in the future. In practice, we expect that only a few 
constraints in the hierarchy need to be examined on average. 

We compute a hierarchical decomposition of the front, say using METIS dKarl |KK98l l as a mesh partitioning 
tool or using a A;d-tree JdBSvKOOOI l. We build the cone hierarchy bottom-up, combining pairs of bounding cones by 
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computing a tight bounding cone for them. Each node in the hierarchy (a binary tree) stores a cone tightly containing 
the cones of its children. 

The cone hierarchy is traversed top-down. Starting at the root, we maintain a priority queue of cone constraints. 
At each step, we query the strictest constraint from the queue. If the current strictest cone is a leaf in the hierarchy 
or if it guarantees sufficient progress (say a constant fraction of the amount of progress when limited by the causality 
constraint alone), then we accept the resulting tentpole height. Otherwise, we replace the strictest cone by its children 
and repeat. 

Refinement and coarsening of the front induce a corresponding refinement and coarsening of the cone hierarchy. 
Since the front adaptivity in response to the error indicator can be be non-uniform, the cone hierarchy can get heavily 
unbalanced. To maintain efficiency, the cone hierarchy needs to be periodically rebalanced or recomputed. 

NONCIRCULAR BOUNDING CONES Anisotropic cone constraints, i.e., cones with noncircular cross sections, com- 
plicate the data structures used for the cone hierarchy but do not introduce substantial algorithmic difficulties. For 
instance, if cones of influence are elliptical, then at each step while constructing the hierarchy we need to compute 
an elliptical cone that contains as tightly as possible the two elliptical cones of its children. Likewise, queries of the 
cone hierarchy are a little different — each requires checking whether a triangle or cZ-simplex of the front intersects 
an elliptical cone. We are able to exploit the facet that a cone of any cross-section is a ruled surface to reduce the 
dimensionality of the simplex involved in checking for intersections. 

Parallel data structures The Tent Pitcher algorithm and all algorithms in this thesis are inherently parallel — 
patches created by pitching non-adjacent local minima of the same front can be solved simultaneously and indepen- 
dently of each other because patch boundaries are causal. 

The front is decomposed into subdomains by a hierarchical partition. Within each subdomain, a local neighbor- 
hood A^ can be advanced independently of other subdomains unless A^ is at the boundary of its subdomain in which 
case interprocessor communication is required to ensure mutual exclusion. The bounding cone hierarchy within the 
subdomain assigned to a single processor is a local sequential data structure. We need a parallel tree data structure to 
query bounding cones stored on other processors. 

Load balancing means repartitioning of the domain. The bounding cone hierarchy must also be recomputed. 

Lazy updates to reduce interprocessor communication At every step, the wavespeed is computed 
as part of the solution and the bounding cone hierarchy needs to be updated to reflect this change. An increase in the 
wavespeed somewhere may impact the progress at a remote location. Therefore, this change propagates potentially to 
every other subdomain, requiring communication and synchronization between processors. 

We can reduce interprocessor communication by performing lazy updates of the bounding cone hierarchy to reflect 
changes in the wavespeed. Because of the no-focusing assumption, current estimates of the wavespeed are always 
conservative and therefore acceptable approximations of future wavespeed. If sufficient progress is guaranteed even 
with these conservative approximations, then there is no immediate need to compute the new wavespeed. Therefore, 
we can queue the updates of the cone hierarchy and perform the actual update later. 
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Chapter summary 

We have been able to use the new operations in our arsenal, such as pitching inclined tentpoles and performing edge 
flips, to give advancing front meshing algorithms in ID- and 2D x Time. Extending the heuristics of this chapter 
to higher dimensions, say 3D x Time would be very interesting. These newer algorithms remain heuristics in the 
sense that we are not able to prove yet that we will be able to successfully mesh the spacetime domain of interest with 
nondegenerate elements for a general class of motion of the boundary or for tracking a fairly diverse class of spacetime 
features. For instance, it remains an open question whether we can align causal and implicit facets of the mesh with an 
arbitrary set of piecewise linear constraints in spacetime. The main challenge remains to prove nontrivial worst-case 
lower bounds on the temporal aspect ratio of spacetime elements under these additional constraints. 
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Chapter 8 

Conclusion and open problems 



We extended Tent Pitcher, the advancing front algorithm due to Ungor and Sheffer JUSOOl) and Erickson et aJ. (IEGSU02| |. 
to adaptively compute efficient, near-size-optimal spacetime meshes suitable for DG solutions that adapt to a posteriori 
spacetime error indicators as well as to nonlinear and anisotropic physics. Our primary motivation was to prove that 
such an algorithm was at all possible with provable theoretical guarantees. In addition to our theoretical results, an 
important contribution of our work is that the algorithms are currently being fruitfully implemented and tested on prob- 
lems of practical complexity. We were able to unify earlier results that considered nonlocal cone constraints ( IThi04l) 
and mesh adaptivity jACE^Oll separately. Our extensions to Tent Pitcher retain the benefits of earlier algorithms such 
as ease of implementation and inherent high degree of parallelism. Additionally, we aim to solve robustness issues 
and generate even better quality meshes in practice than those guaranteed in theory. 

Several theoretical problems remain open and are actively being researched by this author and others. We discuss 
them in the next and final section. 

8.1 Open problems 

There are many exciting avenues to explore and many more problems are likely to be discovered as our algorithm 
are used to solve increasingly realistic problems. In this section, we highlight some promising directions for future 
research. 

Adaptivity in arbitrary dimensions An important open problem is extending adaptive refinement and 
coarsening to spatial dimensions li > 3. Extensions of newest-vertex bisection to higher dimensions are known jAMPOOl 
|B91| ). However, it is challenging to devise necessary and sufficient progress constraints that guarantee positive 
progress. It is also important that the constraints be easy to implement. We consider the problem of incorporating 
adaptive refinement and coarsening in 3D x Time into our meshing algorithm to be the most interesting and practically 
relevant open problem. 

Element quality Consider another measure of element quality: the ratio of inradius to circumradius. The 
circumradius to inradius ratio of an element is defined after scaling the time axis by the wavespeed. A larger inra- 
diusxircumradius ratio means a "better" element. Experiments suggest that the inradius to circumradius ratio is a 
useful measure of element quality. Our algorithms guarantee a lower bound on the worst-case ratio as a result of 
causality and progress constraints, but this bound is not very good. At the same time, some elements with small inra- 
dius to circumradius ratio also seem to have small temporal aspect ratio. It is an important research question to devise 
an algorithm to explicitly maximize the smallest inradius to circumradius ratio. 
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Scale the time axis relative to space by the local wavespeed. Now, the spacetime domain, at least locally within 
each patch is identical to {d+ 1) -dimensional Euclidean space. After this scaling, define the quality of a spacetime 
element as the ratio of its inradius to its circumradius — a larger ratio means a better quality element. 

It is clear that maximizing the tentpole height is not hkely to maximize the worst-case quality of spacetime ele- 
ments. The following is an attempt to describe another choice of location for the top of the tentpole so as to guarantee 
element quality. 

Let T be a front. Let f be a local minimum of T. Let PQR be an arbitrary triangle incident on P. Let p denote the 
circumradius of APQR. Any tetrahedron with PQR as a facet must have circumradius at least p. Attempt to create a 
spacetime tetrahedron over PQR with circumradius not much larger than p but with a guaranteed lower bound on its 
inradius. 

For each triangle PQR, consider its diametral circumball but with radius scaled by (1 + 5) for some 5 > 0. We wish 
to place the top P' of the tentpole inside the intersection of all such circumballs so as not to increase the circumradii 
of resulting tetrahedra by more than (1 + 5). 

In addition, we wish to ensure that the volume of each tetrahedron is bounded from below. We do this by ensuring 
that the distance of P' from the plane of APQR is bounded from below. Thus, we place P' above the plane at an 
orthogonal distance of 7 times the width of APQR. 

Overall, we need to ensure that the intersection of all the circumballs and all halfspaces is nonempty so that we 
have at least one candidate placement of P', the top of the tentpole at P. Note that the final tentpole need not be 
vertical. 

Empirical study of heuristics Even for the basic linear nonadaptive algorithm of Chapter|2] an experimental 
study is needed to try several heuristics for which portion of the front to advance in each step, i.e., which local minimum 
vertex to pitch next when several candidates are available. Different heuristics have been observed to affect the mesh 
quality differently in ID x Time. With the more general advancing front operations of Chapter |7j it is important to 
prioritize various operations that modify the front to ensure that the quality and efficiency of the mesh is improved in 
practice. 

Provably correct boundary tracking We would like to devise a provably coiTect and complete boundary 
tracking algorithm for interesting classes of motion. In many practical situations, the topology of the domain also 
changes with time ( lfPCS99b . For instance, a rocket fuel grain may split into multiple components due to combustion. 
We would be interested in an advancing front algorithm, initially in 2D x Time, that handles changes in the topology of 
the domain over time. In higher dimensions, a wider array of front operations is possible and necessary for boundary 

tracking. For instance, edge flips in 2D x Time generalize to 2 3 flips that replace two adjacent tetrahedra by 

three tetrahedra sharing an edge, and vice versa. The problem of devising a boundary tracking algorithm in arbitrary 
dimensions remains open and of a great deal of interest even in 3D x Time. 

Non-greedy algorithms It should be clear that the progress constraints can be modified in different ways to 
meet different mesh generation objectives. Research so far has focused on worst-case guarantees for temporal aspect 
ratio. Much more research is needed to guarantee worst-case quality in the sense of inradius to circumradius ratio. 
Even more desirable is to guarantee a certain distribution of elements by quality measure. We believe that greedy 
algorithms will not suffice to guarantee worst-case inradius to circumradius ratio in theory and will not perform well, 
without additional heuristics, in practice. 
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Another context in which greedy algorithms are unlikely to fare better than current work is in the case of tracking 
moving boundaries. Intuition suggests that we need an algorithm that constructs only "robustly progressive" fronts. 
We say that a front T is robustly progressive if T is progressive and if T has a local minimum p such that for every 
sufficiently small but finite motion of p with finite velocity the resulting front f is also progressive. For a front T 
to be robustly progressive, it suffices that T has a local minimum p such that (i) for every neighbor q of p, we have 
t(^) — t{p) is bounded away from zero, and (ii) the kernel of link(/?) contains a disk of finite radius centered at p. 
Deriving such an algorithm is an open question. 

NON-SIMPLICIAL ELEMENTS To perform more general front advancing operations, those that change the com- 
binatorics of the spatial projection as well as advance the front in time, it seems necessary to allow more general 
spacetime elements in the near future. Satisfying the two requirements that the spacetime mesh be a weak simplicial 
complex and that every element have at least one causal outflow facet will be harder with the incorporation of more 
general operations. The numerical techniques can be extended to other linear elements like hexahedra, prisms, and 
pyramids with very little difficulty. The data structures used to represent a patch (admittedly, only a transient object) 
would be a little more complicated because incidences between elements and their facets, and adjacencies between 
elements will be a little more complicated. I anticipate no significant difficulties in extending the advancing front 
framework to create and solve patches of mixed types of linear elements. 

More complicated patches Our current algorithms advance a local neighborhood N of the front T to a new 
neighborhood A^' of the new front t'. How complicated can the neighborhoods A^ and A^' be? Can we triangulate the 
spacetime volume between A^ and A^' with a small number of spacetime elements? Given the spacetime volume, can 
we decompose it into a minimum number of good quality spacetime elements? There is a tradeoff between the time 
to solve a single patch versus the number of patches in the mesh for a given spacetime volume. On the one hand, 
we would like patches of small complexity so that the time to solve one patch, even for high polynomial order, is not 
large. This is especially important for an incremental solution strategy where the physical parameters that govern the 
solution are changing rapidly. On the other hand, we would like to minimize the total number of patches for a given 
accuracy to reduce the total computation time. 

The number of elements in a patch created by pitching a vertex p is equal to the number of simplices in star(/?). 
We are already considering operations, such as edge dilation and edge contraction, where the number of elements in 
a patch could be twice as much. The larger is the neighborhood advanced in a single step, the less is the degree of 
parallelism of the meshing algorithm. 
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